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SECTION  1 


INTRODUCTION 


The  probabilities  of  first  aircraft  failure  under  a 
prescribed  mission  loading  and  a number  of  rigorous  and  cursory 
inspections  are  calculated  analytically.  Three  programs  were 
written  for  this  purpose.  The  first  program,  called  Maneuver 
Loading  (MANLOD) , calculates  the  average  of  rise  and  fall  of  the 
stress  process  due  to  aircraft  maneuvers  or  gust  utilizing  a 
simulation  technique  of  random  processes.  The  results  from  this 
program  serve  as  an  input  to  the  other  two  programs.  The  second 
program,  called  Aircraft  Reliability  when  Time  to  Crack  Initia- 
tion is  a Random  Variable  (ARELl) , evaluates  probabilities  of 
first  aircraft  failure  based  on  the  condition  that  time  to  fa- 
tigue crack  initiation  is  distributed  according  to  a two-parameter 
Weibull  distribution.  The  third  program,  called  Aircraft  Rel- 
iability when  Initial  Cracks  are  Random  Variables  (AREL2 ) , 
calculates  probabilities  of  first  aircraft  failure  based  on  the 
initial  crack  size  distribution  model.  The  following  sections 
of  this  document  discuss  these  three  programs,  their  options, 
the  input  data,  and  the  output. 
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SECTION  2 


MANEUVER  LOADING  PROGRAM  (MANLOD) 

The  purpose  of  MANLOD  program  is  to  determine  the  rise 
and  fall  of  the  stress  response  process  due  to  aircraft  maneu- 
vers or  atmospheric  turbulence.  The  results  from  this  program 
serve  as  an  input  to  the  aircraft  reliability  programs  ARELl 
and  AREL2  which  are  described  in  Sec.  3 and  4.  A simplified 
flow  chart  of  the  MANLOD  program  is  shown  in  Figure  1. 

2.1  SIMULATION  AND  RANGE  COUNT 

The  time  history  of  stress  response  due  to  aircraft 

maneuvers  or  gust  is  generated  from  the  specified  stress  spectral 

densities  utilizing  simulation  techniques  of  Gaussian  random 

processes.  For  the  purpose  of  this  study,  the  stress  history 

due  to  gust  is  generated  from  a specified  normalized  spectral 

density  function  corresponding  to  a truncated  Gaussian  white 

noise  spectra.  The  area  under  this  spectral  density  is  unity. 

The  width  of  this  spectra  is  controlled  by  a parameter  called 

BETA  which  ranges  0.0  < BETA  < 1.  For  example,  BETA  = 0 

corresponds  to  a wide  band  process  and  BETA  = 0.9  to  a narrow 

band  process.  The  stress  time  history  due  to  aircraft  maneuvers 

is  generated  from  the  specified  spectral  density  in  Function 

FUNSPT.  If  the  spectral  density  due  to  maneuvers  is  specified 

2 

in  terms  of  accelerations  g , the  conversion  to  stress  spectral 
density  is  achieved  through  the  parameter  g.  To  account  for 
the  unsymmetric  exceedance  observed  in  actual  maneuver  exceedance 
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START 


END 


Figure  1 


Simplified  Flow  Chart  of  MANLOD  Program 


curves,  the  stress  time  history  is  modified  by  multiplying  the 
negative  portion  by  a factor  called  REDUCE  = standard  deviation 
for  negative  maneuvers/standard  deviation  for  positive  maneuvers. 

The  rise  and  fall,  As^,  of  the  simulated  process  S (t)  are 
determined  using  two  separate  range  count  methods.  In  the  first 
method,  the  AS-,  . values  are  determined  by  selecting  the  difference 

-L  # 1 


second  counting  method,  the  AS„  . values  are  selected  from  the 

Z r ± 

difference  between  those  peaks  and  valleys  which  are  separated 
by  the  Ig  line  as  shown  below. 
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2 . 2 INPUT  DATA 


There  are  no  input  data  cards  used  for  the  Maneuver 
Loading  Program.  The  necessary  data  parameters  need  to  be  spe- 
cified at  the  beginning  of  the  main  program.  The  input  spectral 

2 

density  function  of  the  maneuver  loading  in  g units  need  to 
be  specified.  An  option  is  provided  to  use  a normalized  trun- 
cated white  noise  spectrum.  The  input  parameters  are  as  follow: 
maneuver  loading  spectral  density  function  SS (W)  with  parameters 
A,  FB  and  PHIO  where  A = a shape  parameter;  FB  = break  frequency 
in  Hz;  PHIO  = peak  value  of  the  spectral  density  function  SS (W) ; 
IX  = an  odd  integer  with  less  than  9 digits  where  by  changing 
IX  a new  sample  function  is  generated;  N = number  of  frequency 
subdivisions  on  the  spectral  density  function  SS (W)  in  powers 
of  2;  NPT  = a number  of  simulated  points  of  the  output  sample 
function  specified  by  NPT  = 2 **  Ml  with  2 * N being  equal  to 
or  less  than  NPT,  WU  = the  upper  cut-off  frequency  on  SS (W) ; 

WL  = the  lower  cut-off  frequency  on  SS (W) ; Ml  = a positive  inte- 
ger such  that  NPT  = 2 **  Ml,  NSIGN  = an  option  parameter  where 
for  NSIGN  = 0 a normalized  white  noise  spectral  density  with 
a unit  standard  deviation  is  provided  by  the  program,  and  for 
NSIGN  = 1 the  user  must  supply  a spectral  density  for  the  man- 
euver loading  or  use  the  one  specified  in  the  subroutine  FUNSPT; 
REDUCE  = ratio  of  the  standard  deviations  such  that  REDUCE  = 

a /a  • where  o is  the  standard  deviation  determined  from  the 
n p n 

time  history  below  the  lg  level  and  a ^ is  the  standard  deviation 
obtained  from  the  time  history  above  the  lg  level;  NP  = number 
of  points  plotted  by  the  subroutine  PLOT;  BETA  = parameter  which 
specifies  the  width  of  the  truncated  Gaussian  white  noise  spec- 
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trum;  NSAMPL  = number  of  samples  to  be  simulated. 


2.3  PROGRAM  DESCRIPTION 

The  Maneuver  Loading  Program  consists  of  a main  program 
unit  and  fourteen  subroutines.  The  program  listing  is  given 
in  Appendix  A.  The  main  program  unit  called  MANLOD  and  the 
fourteen  subroutines:  AREA,  SIMULT,  PLOT,  PICK,  STDEVI,  MEAN7 , 

PICKM,  PICKMX,  CROSS,  SIMPSN,  FUNSPT,  HARM,  RANDU,  STATIS  are 
described  as  follows. 

2.3.1  The  Main  Program  Unit  MANLOD 

The  main  program  initializes  the  input  data  and  with  a 
call  to  the  subroutine  AREA  calculates  the  area  under  the  man- 
euver loading  spectral  density  specified  in  Function  FUNSPT. 

The  time  realizations  of  the  loading  (or  stress)  are  generated 
by  calling  subroutine  SIMULT.  When  NSIGN  = 0,  stress  time  his- 
tory corresponding  to  gust  input  is  provided.  When  NSIGN  = 1, 
stress  time  history  due  to  maneuvers  is  given.  The  rise  and 
fall  of  the  simulated  time  history  are  determined  by  calling 
subroutine  PICK  for  the  first  counting  method  and  subroutine 
PICKM  for  the  second  counting  method.  The  rise  and  fall  values 
are  raised  to  an  arbitrary  power  B ranging  from  1 to  10  by  sub- 
routine STDEVI.  Then,  the  average  of  rise  and  fall  raised  to 
B power  is  calculated  by  calling  subroutine  MEAN 7 . The  total 

number  of  peaks  of  the  simulated  process  S (t)  is  counted  by  the 

m 

subroutine  PICKMX  and  the  crossing  rate  is  determined  from  sub- 
routine CROSS.  A subroutine  PLOT  is  called  to  plot  the  output 
time  history  for  one  realization. 
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2.3.2  Subroutine  AREA  (S,  N,  WUP,  WL) 


Subroutine  AREA  calculates  the  area  under  the  curve 
specified  in  FUNSPT.  It  utilizes  a Simpson  rule  type  integra- 
tion procedure  given  in  Function  SIMPSN.  The  arguments  for  this 
subroutine  are 

S = input  data  from  FUNSPT; 

N = number  of  subdivisions  on  spectral  density  SS; 

WUP  = upper  cut-off  frequency  on  SS; 

WL  = lower  cut-off  frequency  on  SS. 

2.3.3  Subroutine  SIMULT  (A,  N,  NPT,  WU,  WL, 

IX,  Ml,  NSIGN,  INV,  S) 

The  subroutine  SIMULT  is  used  to  generate  samples  of  the 
time  history  corresponding  to  gust  or  maneuver  loading.  The 
program  utilizes  the  Fast  Fourier  Transform  (FFT)  algorithm 
by  calling  subroutine  HARM.  The  essential  equations  for  this 
simulation  are 

Real  A ( 2 * J + 1)  = SQRT(2.0  * SP)  * COS  ( PHI) 

Imaginary  A(2  * J + 2)  = SQRT(2.0  * SP)  * SIN (PHI) 

where  J = 1,  2,  . . . , N,  SP  = specified  spectral  density, 

PHI  = random  phase  angles  uniformly  distributed  between  0 and 

2ir  which  are  generated  by  subroutine  RANDU.  The  arguments  for 

this  subroutine  are  A = specified  by  the  above  equations;  N = 
number  of  subdivisions  on  spectral  density  SP;  NPT  = number 
of  time  points  to  be  simulated;  WU  = upper  cut-off  frequency  on 
SP;  WL  = lower  cut-off  frequency  on  SP;  IX  = a number  specified 
in  the  input;  NSIGN  = specified  in  the  input;  INV,  S = internal 
parameters  of  the  subroutine  . 
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2.3.4  Subroutine  HARM  (A,  M,  INV,  S,  IFSET , IFERR) 


Subroutine  HARM  perform  a Fast  Fourier  Transform  (FFT) 
on  the  specified  function  A(I),  (I  = 1,  2,  . . . , N) . The  argu- 

ments for  this  subroutine  are 

A = generated  input  by  the  subroutine  SIMULT; 

M = number  of  points  to  be  simulated,  i.e.,  for  one 
dimensional  case  M(l)  = Ml,  M(2)  = 0,  M(3)  = 0. 

INV,  S = internal  parameters  of  the  subroutine; 

IFSET  = 1,  implies  a forward  Fourier  Transform,  IFSET  = -1 
implies  an  inverse  Fourier  Transform; 

IFERR  = internal  parameter  of  the  subroutine. 

2.3.5  Subroutine  PLOT  (NP,  DT,  A,  0.0,  M) 

Subroutine  PLOT  plots  a given  function  A with  a dimen- 
sion NP  at  step  intervals  DT.  The  argument  M indicates  how  many 
points  from  given  NP  needs  to  be  plotted.  For  example,  if 
M = 1,  it  implies  that  every  point  is  plotted,  if  M = 2,  it  implies 
that  every  second  point  is  plotted. 

2.3.6  Subroutine  PICK  (A,  NPT , S,  NPK(J)) 

Subroutine  PICK  determines  the  rise  and  fall  AS  cor- 

li 

responding  to  the  first  counting  procedure  as  described  in  Sec. 

2.1.  The  arguments  for  PICK  program  are 

A = input  array  from  which  rise  and  fall  needs  to  be 
found ; 

NPT  = number  of  points  in  the  input  array  A; 

S = an  array  in  which  the  rise  and  fall  are  stored; 

NPK  = number  of  rise  and  fall  values  found  from  the  specified 
input  A. 
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2.3.7  Subroutine  STDEVI  (NST , NPK(J),  S,  STD, 

RMS,  A1(J),  A2(J),  . . .,  A5(J),  AN) 

This  subroutine  raises  the  rise  and  fall  values  calculated 
in  PICK  to  powers  Al,  A2,  . . .,  A5,  where 

NST  = parameter  calculated  by  the  main  program; 

NPK  = number  of  rise  and  fall  values; 

S = rise  and  fall; 

STD,  RMS,  AN  = parameters  calculated  inside  the  subroutine. 

2.3.8  Subroutine  MEAN 7 (NSAMPL,  Al,  . . .,  A5, 

AMI,  . ♦ . , AM5 , AA,  AN) 

Subroutine  MEAN 7 calculates  the  average  of  rise  and  fall 
raised  to  arbitrary  power  in  subroutine  STDEVI. 

NSAMPL  = sample  size  for  which  average  is  calculated; 

Al,  . . . , A5  = parameters  obtained  in  subroutine  STDEVI; 

AMI,  . . , AM 5 = calculated  in  the  subroutine; 

AA  = calculated  in  the  main  program; 

AN  = calculated  in  the  subroutine 

2.3.9  Subroutine  PICKM  (A,  NPT,  S,  MPK(J)) 

This  subroutine  calculates  the  rise  and  fall  AS  corres- 

2i 

ponding  to  the  second  counting  procedure  as  described  in  Sec.  2.1. 
The  meaning  of  the  arguments  is  the  same  as  for  subroutine  PICK. 

2.3.10  Subroutine  PICKMX  (Y,  NPT,  PMAX,  LPEAK , 

PMIN , MTROG) 

Subroutine  PICKMX  picks  all  the  peaks  and  troughs  from 
a given  process. 

NPT  = number  of  points; 

Y = input  process; 
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PMAX  = array  for  peak  values; 

PMIN  = array  for  trough  values; 

LPEAK  = number  of  peaks  found; 

MTROG  = troughs  found. 

2.3.11  Subroutine  CROSS  (STDEV,  PMAX, 

PMIN,  LLL , MMM,  VZERO,  RATIO) 

This  subroutine  determines  the  crossing  of  the  threshold 
level  specified  by  STDEV  for  a given  process. 

STDEV  = standard  deviation  of  the  given  process. 

PMAX  = peak  value  array  of  the  process; 

PMIN  = trough  value  array  of  the  process; 

LLL  = number  of  peaks; 

MMM  = number  of  troughs; 

VZERO  = indicates  zero  threshold  crossing; 

RATIO  = when  RATIO  = 1,  crossing  for  the  total  given  time 

history  is  determined,  when  RATIO  = 0,  crossing  per 
hour  is  found. 

2.3.12  Subroutine  SIMPSN  (AF,  NPT , 

DSTEP , SIMP) 

Subroutine  SIMPSN  integrates  a given  function  AF  utilizing 
a Simpson  rule  integration  procedure. 

NPT  = number  of  points  to  be  integrated  (must  be  an  odd 
number) 

DSTEP  = integration  step; 

SIMP  = resulting  area  of  integration. 

2.3.13  Subroutine  RANDU  (IX,  IY,  YFL) 

Subroutine  RANDU  generates  independent  random  numbers 


-10- 


uniformly  distributed  between  0 and  1.  By  multiplying  these  num- 
bers by  2tt,  random  phase  angles,  which  are  necessary  for  the 
digital  simulation  procedure,  are  obtained.  An  odd  integer,  IX, 
needs  to  be  specified  in  the  main  program.  By  changing  IX,  a 
new  seed  of  random  numbers  is  generated. 

IX  = any  odd  integer  less  than  9 digits; 

IY  = internal  parameter  in  the  program; 

YFL  = output  - random  numbers  between  0 and  1. 

2.3.14  Subroutine  FUNSPT  (SP,  Wl) 

Subroutine  FUNSPT  is  used  to  array  the  normalized  empir- 
ical spectral  density  function  corresponding  to  aircraft  maneuvers. 
SS  (W)  = 2tt  * PHIO/  ( 1 + (2tt  * W/FB)  **  A)  /ANORML 

where 

W = frequency  in  HZ; 

FB,  A,  PHIO  = data  parameters; 

ANORML  = area  under  SS  which  is  determined  in  subroutine 
AREA. 

2.3.15  Subroutine  STATIS  (A,  NP,  STD,  RMS) 

Subroutine  STATIS  determines  the  standard  deviation  and 
rms  of  maneuver  loading. 

A = input  array; 

NP  = number  of  points; 

STD  = output  - standard  deviation; 

RMS  = output  - root  mean  square. 

2 . 4 PROGRAM  OUTPUT 

The  MANLOD  program  generates  printed  output.  In  Table 
1/  the  average  of  rise  and  fall  for  one  realization  is  given 
using  the  first  range  counting  method  described  in  Sec.  2.1. 
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TABLE  1,  RISE  AND  FALL 
(FIRST  COUNTING  METHOD) 


3 FT  A = 0.0  F A K = ??85  4J  = 0.2000  V.  = 0.0070 

F A vj  1 = 0.60?  M-AN7  = 0.9^  V'-  S'l’.r  l.SAT  MFAVT5“=  ?.  1 63  ME  A N4  = 3.129 

*1EAV3=  0. 1 7 3 D B F 0?  MFA\|7=  O.43034E  0?  M E A N 8 = 3.13B26E  07  ME  A M 9=  0.42070E  03 

MfAN4fr=  4.661  MFAN5=  7.126 
MEAN  1 0=  0.13425E  04 

STOEV=  0. 6505576E  00  RMS=  0.9499773E  00 

NPEAK  = number  of  rise  and  fall; 

WU  = upper  cut-off  frequency  on  spectral  density; 

WL  = lower  cut-off  frequency  on  spectral  density^- 
STDEV  = standard  deviation  of  the  rise  and  fall; 

RMS  = root  mean  square  of  the  rise  and  fall; 

MEAN1,  MEAN 2 . . . MEAN10  = indicates  the  mean  values  if  rise 

and  fall  which  were  raised  to  powers  1,  2,  . . 10^ 

respectively. 
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Similar  results  are  shown  in  Table  2 for  the  second  range  count- 
ing method.  The  standard  deviation  and  rms  of  S (t)  are  also 

m 

included  in  this  table. 


TABLE  2.  RISE  AND  FALL 
(SECOND  COUNTING  METHOD) 


3 ET  A = 0.0 


E fl.  N 1 = 1.  29? 

'll  AN  3=  3 .3  44}  3 E 


NPEAK=  1000  tfj=  0.2000 

MEA  N 2 = 2.274  ME  AN3=  4.791 

METTT=  T773TTTTE  03 M~A  N3  = 


WL  = 0.0070 


MEAM35- 
7 .'8T775  £~ 


7.  2 99 
T*r~ 


STDEV=  0. 7776743E  00 

MEAN4=  11.401  MEAN45= 
*FAT9=  ~ Ot  29-05-M-  tT4 ME  AMI  0= 


RMS=  0.  15080346  01 

18.237  ME  AN5=  29.802 

05 — — 


TTANBARD  171V1  A TIflN  UP  HE  TIME  HJSHJMV 

? TOE/  = 0 .5485  804 F 00  R . M . S .=  0.  7 13095 8E  00 


The  average  of  rise  and  fall  corresponding  to  averaged  eight 
realizations  is  given  in  Table  3 for  the  first  range  counting 
method  and  in  Table  4 for  the  second  range  counting  method. 
The  number  of  rise  and  fall  for  each  method  is  also  indicated 
in  these  tables. 
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TABLE  3.  AVERAGE  OF  POWER  OF  RISE  AND  FALL 
(FIRST  COUNTING  METHOD) 


CASE  I 


A7E9.A0F  OF  P9WE3  OF  RISE  AMD 

FALL 

Avn  n’j^p.fr  of  rise  ano  hall 

POVE  0 

RFTA 

1.0 

9.0 

775  375  4775  47? 

— 770 

91  STE  F,  'FALL- 

n .0 

_ n . 6 V 6 

0. 933 

1 . 679  ?.  397  3.568  5.4 

79  3.645 

2269 

R R T A 

6.0 

77 

N 870  9.  u 

LU. 

o_.o 

?3.06 

6 6.  90  203.68  659.  41 

2235.61 

2269 

TABLE 

4.  AVERAGE  OF  POWER  OF  RISE 

AND  FALL 

(SECOND  COUNTING  METHOD) 


CAS-  II 


AV^IAS^  n=  POWER  OF  RISE  AMO  FALL 

AMO  NJM3ER  OF  RISE  AMD  pA__ 
P0W2  3 

BETA 

1.0 

2 . 0 

77t  775  470  475  570 

< 1st  6 FALL 

_ Q.0- 

U 7J-6. 

9 . 7 64 

4.371  7.51?  1 1 .932  1 9.3  9 5 32.211 

1031 

- — — ~ r . o 971  TT7  * E A K S' 

3 FT4 

0.  o 9'±,  00 ? 3 2 . 9 6 6 .31 33  50.7?  1 2147 . 2 7 1031 


The  level  of  threshold  crossing  per  hour  against  standard 
deviations  of  the  process  is  given  in  Table  5. 
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TABLE  5.  THRESHOLD  CROSSING 


SIMULATE! 

TOTAL  5 2 OS  SI  'l  GS  AT 
CROSSING  W tr  H BOTH 

0 = S UL  T 

OirpERENT  LEVEL  OF  SIGMA, 
fVF  AND  - V;  SLJR 

- 5. OSTGMA= 

0 

",  R05  5T  MG  5 

-4.8S I GM  A = 

0 

CROSS! NOS 

-4.65  Tf,MA  = 

3 

T,  R0S  5 I MGS 

-4.4SIGMA= 

0 

CROSSI  NGS 

-4.2S  IGMA= 

0 

CROSSINGS 

-4.  0$IGMA= 

0 

CROSS! NGS 

-3.85  TOM A = 

0 

C ROSSI  NGS 

-3.55  TOM  A= 

3 

: R05  5 I NGS 

- 3,4SI3MA= 

0 

: ROSSI  MGS 

-3.25  I GMA= 

0 

GR0SSI NGS 

— 1 1 Wil  Hi  111 

3 

r ROS  5 T NGS 

-2.8S10MA* 

0 

; ROSSI NGS 

-2.6SIGMA= 

0 

CROSSINGS 

-2. 45  T GMA= 

0 

C ROSSI NGS 

3 

kUiLiIMn 

-2.0STGMA= 

0 

C ROS  S I NGS 

-1.8SISMA* 

0 

C ROSSI  NGS 

-1 .65 IGMA= 

0 

CROSSINGS 

-1  . 45  T GM  A = 

. 3 

C ROSS  I NGS 

- 1 »2SI SM A = 

0 

CROSS!  NGS 

-1.3S  TOMA= 

0 

CROSSINGS 

-0.  851  GMA= 

4 

CROSS!  NGS 

wiiwniii 

25 

— Ill  1 

-0.4S I0MA= 

84 

C ROSS  T NGS 

- 0.2 SI 3MA  = 

188 

CROSSI NGS 

0.35 1GMA= 

3 50 

CROSS! NGS 

0.2  SIGMA* 

CROSS  T NGS 

0 . 4S  T GM  A = 

267 

: . : 

0.65 IGM  A= 

229 

0.  85  I GM  A* 

1 93 

CROSSI NGS 

1 .05  TGM A* 

1 63 

C ROSS!  NGS 

1 . 2 S I GM  4 = 

134 

C RoS  5 HOG'S 

l . 4S  T GM  A = 

1 08 

: ROSS!  NGS 

1 . 6 5 I G M A= 

78 

CROSSINGS 

84 

CROSSI NGS 

2. OS IGMA= 

40 

C ROSSI  NGS 

2 . 2 5 T GM  A = 

24 

CROSS  I NGS 

2. 4SI GM A= 

13 

CROSSI  NGS 

8 

?.8SIGMA= 

3 

3.051 GM  A= 

2 

: ROS  ST  NGS 

3.2  5 I GM  A= 

1 

CROSS  I NGS 

3.45T  G M A = 

0 

“ROSSI  NGS 

3.65 IGM  A* 

0 

CROSSI NGS 

3.3SIGMA* 

0 

C ROSS  I NGS 

4. 05 I GM A= 

0 

CROSSI NGS 

i(Ciu«*T~P 

0 

CROSSINGS 

r R.  OS  S I MGS 

CROSSI  MGS 

CROSS  T NGS 

CROSSI MGS 
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Section  3 


AIRCRAFT  RELIABILITY  WHEN  TIME  TO  CRACK  INITIATION 
IS  A RANDOM  VARIABLE  (AREL1) 

The  purpose  of  this  program  is  to  calculate  the  probabil- 
ity of  first  aircraft  failure  either  of  a single  airplane  or  of 
an  airplane  in  a fleet.  The  analytical  probability  of  failure 
model  is  constructed  for  the  case  when  time  to  fatigue  crack  ini- 
tiation at  a critical  location  on  the  aircraft  is  a random  variable 
having  a Weibull  type  distribution.  The  analytical  model  allows 
the  inclusion  of  rigorous  base  type  inspections  which  are  per- 
formed at  equal  time  intervals  over  the  design  life  of  the  air- 
craft. A simplified  flow  chart  of  the  program  is  shown  in  Figure 
2*  The  following  paragraphs  discuss  the  input  data,  the  struc- 
ture of  the  program,  and  the  output.  The  section  on  the  structure 
of  the  program  includes  only  the  key  analytical  expressions. 

For  a detailed  description  of  all  the  equations,  the  reader  should 
consult  the  Technical  Report. 
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START 


INPUT  DATA 


SENSITIVITY  STUDY 


Q 


A (T) 


R(T) 


H (T) 


TH  (T) 


SUBROUTINES 


4 HCOM  (H,  ARO,  AR1 , , AR8 , P.  AN) 

HSIN  (H,  ARO,  AR1 , AR8 , P.  AN.  UPC) 


1 


I 


L 


QSF  (H,  Y,  Z,  NDlM) 


PROBABILITY  OF  FAILURE 
NO  INSPECTION 


PROBABILITY  OF  FAILURE 
WITH  INSPECTION 
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r 


ERF (X) 

-I EEXP  (X) 


F (X) 
■FS (X) 


END 


Figure  2.  Simplicied  Flow  Chart  of  the  AREL1  Program 
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3.1  INPUT  DATA 


There  is  one  input  data  card  used  for  the  AREL1  program. 
The  Format  of  this  data  card  is  illustrated  in  Table  6. 

TABLE  6. 

DATA  CARD  FORMAT  FOR  AREL1  PROGRAM 


Column 

Contents 

Format 

1-10 

XY  = a number  from  1 to  32 
which  references  the  32 
parameters  in  the  system 

F 10.0 

11-20 

‘ 

YZ1  = new  input  data  if 
only  one  parameter  is 
changed 

F 10.0 

21-30 

YZ2  = new  input  data  if  two 
parameters  are  changed 

F 10.0 

31-40 

YZ3  = new  input  data  if  three 
parameters  are  changed 

F 10.0 

The  remaining  data  parameters  are  specified  at  the 
beginning  of  the  main  program.  For  the  physical  description 
of  these  parameters  see  the  Technical  Report. 


The  data  are  as  follow:  EF  = A , material  constant;  BE  = b, 

a parameter  in  the  fatigue  crack  propagation  model;  DKTH  = 
threshold  value  of  the  range  of  stress  intensity  factor; 

FED  = X',  material  constant,  C = c,  material  constant:  KC  = 
critical  stress  intensity  factor;  AO  = initial  fatigue  crack 
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size;  ALPHA  = shape  parameter  in  Weibull  distribution;  BETA  = 
scale  parameter  in  Weibull  distribution;  CMP(l)  = 1.,  CMP(2)  = 

1.,  CMP (3)  = 1.  implies  that  the  aircraft  loads  due  to  nonstorm 
turbulence,  thunderstorms,  and  maneuvers  each  is  taken  as  com- 
posite Gaussian  processes;  CMP(l)  = 0.0,  CMP (2)  = 0.0,  CMP (3)  = 
0.0  implies  that  the  aircraft  loads  due  to  nonstorm  turbulence, 
thunderstorms,  and  maneuvers  each  is  taken  as  a single  Gaussian 
process;  UPC(l)  = 1.,  UPC (2)  = 1.,  UPC (3)  = 1.  implies  that  the 
input  loads  are  truncated  at  certain  level  (single  Gaussian 
case  only);  EL  = the  truncation  index  on  maneuver  loads; 

G = gravitational  constant  in  ksi;  AA(4),  AAA(4)  = rise  and 
fall  parameters  of  the  random  stress  process  due  to  gust  which 
can  be  determined  from  MANLOD  program  by  selecting  the  option 
NSIGN  = 0;  AAAA ( 4 ) = rise  and  fall  parameter  of  the  random 
stress  process  due  to  maneuvers  which  can  be  determined  from 
MANLOD  program  by  selecting  the  option  NSIGN  = 1;  BAB (4)  = 
(2*1-1)!!  where  I = BE/2  and  (2  * I - 1) ! ! = 1 • 3 • 5 . . . ; 
AN(1),  AN(2)  = characteristic  stress  response  frequencies 
corresponding  to  gust;  AN (3)  = characteristics  stress  response 
frequency  corresponding  to  maneuvers;  ANA  = characteristic  stress 
response  frequency  corresponding  to  ground  - air  - ground 
loads;  PP(1),  PP(2),  PP(3)  = fractions  of  aircraft  flight  time 
spent  in  clear  air  turbulence,  thunderstorms,  and  maneuvers, 
respectively;  SGC(l),  SGC(2),  SGC(3)  = stress  intensities  due 
to  clear  air  turbulence,  thunderstorms,  and  maneuvers,  respec- 
tively, when  each  loading  is  treated  as  a composite  Gaussian 
process;  SGS(l),  SGS(2),  SGS(3)  = stress  intensities  due  to 
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clear  air  turbulence,  thunderstorms,  and  maneuvers,  respective- 
ly, when  each  loading  is  treated  as  a single  Gaussian  process; 
XO  = 1G  level  loading  (level  flight) ; ZR  = stress  intensity 
due  to  ground-air-ground  loads;  VO  = coefficient  of  mater- 

ial strength  variation;  AMVO  and  RO  = mean  of  the  material 
strength;  FSF  =1.0  indicates  a fail-safe  type  design,  FSF  = 

0.0  indicates  a slow  crack  propagation  model;  ZETA  = factor 
indicating  the  residual  strength  at  the  fail-safe  crack  size 
AS;  AS  = fail-safe  crack  size;  TT  = design  service  life  of  the 
airplane;  A1  = minimum  crack  size  below  which  the  crack  cannot 
be  detected  by  inspection;  A2  = maximum  crack  size  beyond  which 
the  crack  is  always  detected  by  inspection;  RM  = power  expon- 
ent for  the  crack  detection  probability  model;  U = inspection 
probability;  NN  = number  of  rigorous  base  type  inspections; 

NPT  = number  of  divisions  on  Weibull  distribution;  NOI  = number 
of  output  points  when  no  inspections  are  included  in  the  air- 
craft probability  of  failure  model;  NF  = number  of  aircraft 
in  a fleet.  The  units  of  the  input  data  are  time  = hours, 
length  = inches,  strength  = ksi;  force  = pounds. 

3.2  PROGRAM  DESCRIPTION 

The  Aircraft  Reliability  Program  ARELl  consists  of  a 
main  program  unit  and  seven  subroutines  HCOM,  HSIN,  QSF,  ERF, 
EEXP , F and  FS  which  are  described  in  the  following  paragraphs. 
The  program  listing  is  given  in  Appendix  B. 

3.2.1  The  Main  Program  Unit  ARELl 

The  main  program  reads  and  initializes  the  specified 
data  as  described  in  Sec.  3.1.  The  first  part  of  the  program 
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is  designated  for  the  sensitivity  study  of  the  32  parameters  in 
the  aircraft  reliability  model.  This  is  achieved  by  specifying 
the  value  of  XY  from  1 to  32.  Each  number  of  XY  corresponds 
to  a particular  parameter  in  the  system  which  are  defined  in 
Sec.  3.1. 

The  input  load  parameter,  Q,  due  to  the  combined  effect 
of  the  clear  air  turbulence , thunderstorms,  maneuvers,  and 
ground  - air  - ground  loads  is  calculated  from 

A = AN ( 1 ) * PP(1)  * SB  (1)  + AN ( 2 ) * PP(2)  * SB(2) 

+ AN (3)  * PP(3)  * SB ( 3 ) + ANA  * ZBB 

where  SB(1),  SB(2),  SB(3)  and  ZBB  are  the  averages  of  the  BE  power 
of  the  rise  and  fall  of  the  stress  response  process  corresponding 
to  clear  air  turbulence,  thunderstorms,  maneuvers,  and  ground  - 
air  - ground  loads,  respectively.  These  parameters  are  evaluated 
for  a composite  Gaussian  process  and  for  a single  Gaussian 
process.  The  fatigue  crack  size  at  flight  time  T,  A(I) , (I  = 

1,  2,  . . .,  NPT)  is  determined  utilizing  the  fail-safe  design 

and  the  slow  crack  propagation  models.  The  residual  strength, 

R ( I ) , is  determined  from  the  Griff ith-Irwin  equation 
R ( I ) = FLOAT (KC)  * SQRT ( 2 . 0/PI/A ( I ) ) /RO 
for  the  slow  crack  propagation  model  and  from 

R ( I ) = 1.0  - (1.0  - ZETA)  * SQRT ( (ABS (A (I) 

- AO)  / (AS  - AO)  ) ) 

for  the  fail-safe  design,  where  RO  = ultimate  strength.  Then, 
the  failure  rate,  H(I),  is  determined  by  calling  subroutine 
HCOM  when  the  input  loading  is  treated  as  a composite  Gaussian 
process  and  subroutine  HSIN  when  the  input  loading  is  taken  as  a 
single  Gaussian  process.  The  summation  of  the  failure  rate  from  the 
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fatigue  crack  initiation  to  T flight  hours  after  crack  initiation 
is  obtained  by  integrating  the  failure  rate  H (I)  by  calling  sub- 
routine QSF . The  time  to  crack  initiation  having  a two  para- 
meter Weibull  distribution  is  calculated  from 

W(I)  = ALPHA/BETA  * TBETA  **  ALPHl 
* EXP ( -TBETA  **  ALPHA) 

where  TBETA  = T/BETA,  ALPHl  = ALPHA  - 1.  Then,  by  specifying 
NN  = 0 , which  implies  no  inspections,  the  probability  of 
first  aircraft  failure,  P(I),  is  calculated.  When  NN  / 0 with 
the  inspection  probability  F specified  in  FUNCTION  F,  the  pro- 
bability of  first  aircraft  failure,  P(J),  (J  = TT/NN)  is  deter- 
mined at  each  inspection  interval.  The  probability  of  first 
aircraft  failure,  TP(J),  in  a fleet  of  NF  airplanes  is  calcu- 
lated from 

TP  (J)  = 1.0  - (1.0  - P(J))  **  NF 

3.2.2  Subroutine  HCOM  (H,  ARO , ARl , . . .,  AR8 , P,  AN) 

The  subroutine  HCOM  determines  the  failure  rate  H(I) 
for  the  case  when  input  loading  is  taken  as  a composite  Gaus- 
sian process.  The  necessary  arguments  of  the  program  are 
H = failure  rate  - output; 

ARO,  ARl,  . . .,  AR8  - calculated  in  the  main  program; 

P = input  data  corresponding  to  PP(1),  PP(2)  and  PP(3); 

AN  = characteristic  frequency  calculated  in  the  main  program. 

3.2.3  Subroutine  HSIN  (H,  ARO,  ARl,  ■ . AR8 , P,  AN,  UPC) 

The  subroutine  HSIN  calculates  the  failure  rate,  H(I), 
for  the  case  when  input  loading  is  considered  as  a single  Gaus- 
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sian  process.  The  subroutine  provides  an  option  for  trunca- 
ting the  failure  rate  at  certain  truncation  level  EL  (for 
example  maneuvers).  This  is  achieved  by  setting  UPC  = 1.0 
in  the  main  program.  The  necessary  arguments  of  the  program 
are  the  same  as  in  subroutine  HCOM  and  UPC  is  specified  as 
input  data. 


3.2.4  Subroutine  QSF  (H,  Y,  Z,  NDIM) 

Subroutine  QSF  integrates  the  given  input  function  Y 
at  specified  integration  steps  H and  stores  the  results  in  Z 
at  each  integration  step.  NDIM  = number  of  input  points  to 
be  integrated. 


3.2.5  Function  EEXP  (ARG) 

Function  EEXP  (ARG)  is  used  to  calculate  EEXP  (ARG)  = 
EXP  (ARG)  - 1.0  by  utilizing  an  expansion  procedure  for  the 
case  when  ARG  > -50.0  and  0 < ARG  < .01.  When  ARG  < -50.0, 
EEXP  (ARG)  = -1.0,  and  when  ARG  > .01  no  expansion  is  used. 

3.2.6  Function  F(A) 

Function  F(ft)  specifies  the  crack  detection  probabi- 
lity corresponding  to  rigorous  type  inspections  in  terms  of 
the  crack  size  A.  The  necessary  data  for  this  function  is 
specified  in  the  main  program. 

3.2.7  Function  FS (A) 

Function  FS (A)  determines  probabilities  of  not  de- 
tecting a fatigue  crack. 

3.2.8  Function  ERF 

Function  ERF  calculates  the  error  function  in  terms 
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of  exponentials. 


3 . 3 PROGRAM  OUTPUT 

The  AREL1  Program  generates  printed  output.  First  the 
input  data  XY,  YZ1,  YZ2 , YZ3  is  printed  as  shown  below  for  a 
particular  example- 


XY,YZ  1 ,YZ2,YZ3  0.  3000000F  02  0.  0 0 .3 


0.0 


The  calculated  parameters  SBl,  SB2 , SB3,  ZBB,  Q,  ANS,  ATH,  TB, 


AC,  TC,  CDS,  FSF  are  printed  as  shown  below 
SBl  » SB2  » S B3  t Z BB?  0,  ANS 

0.8283E  02(  1.0)  0.3622E  04(  1.0)  0.2400E  33(  1.0) 


0.3063E  05 


0 .6267E  05  0.3300F  03 


F04  (T.  I 0T51:  00  0 .482  ^ E 04""  0 .1  434b -0  l U.  1U0DE  "0T 

where 

SBl,  SB2,  SB3 , ZBB  - average  stress  intensities  raised  to  BE 
power  corresponding  to  clear  air  turbulence;  thunderstorms, 
maneuvers,  and  ground-air-ground  loads,  respectively. 

Q = combined  input  loading; 

ANS  = PP(1)  * AN (1)  + PP ( 2 ) * AN (2)  + PP(3)  * AN(3)  is  tne 
average  characteristic  frequency  of  gust  and  maneuver  loads 
combined; 

ATH  = fatigue  crack  size  at  the  threshold  level; 

TB  = time  to  reach  threshold  crack  size; 

AC  = initial  critical  crack  size; 

TC  = time  to  reach  AC; 

CDS  = c ; calculated  material  constant  in  the  crack  propagation 


model 


FSF  = index  denoting  fail-safe  or  slow  crack  propagation  model. 


Table  6 
A,  R,  H,  TH  and  ' 


contains  a listing  of  the  output  sample  for 
as  functions  of  flight  time  T. 
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TABLE  6.  CRACK  SIZE  A,  RESIDUAL  STRENGTH  R,  FAILURE  RATE  H 


SUMMATION  OF  FAILURE  RATE  TH,  TIME  TO  CRACK  INITI- 
ATION W. 


I »T  IMEt A »R  » H,TH,W 


1 
51 
101 
151 
201 
251 
301 
351 
401 
451 
5 01 
5 51 
601 
651 
701 
751 
8 01 
851 
901 
9 51 
1 001 
1051 
1 101 
1151 
1201 
1251 
13  01 
1351 
1401 
1451 
1501 
NN= 


O'.  0 
500.0 

1 000.0 

1500.0 

2000.0 

2500.  0 

3000. 0 

3500.0 

4000.0 

4500.0 

5000.0 

5500.0 

6000.0 

6500.0 

7000.0 

7500.0 
"8000.  0 

8500.0 
°000.0 
9500.  0 

TOO 00.0 
1 0500.0 
11000.0 

11500. 0 
12  00  C. 0 

12500.0 

13000.0 

13500. 0 
1400  0 .0 

14500.0 

15000.0 
0 


0. 0400 

0 .0412 
0 .0427 
0. 0445 
0.0467 
0.0496 
0.0536 
0.0603 
C . 0736 
0 .2145 

1 .5828 
7. 0090 
7.0000 
7.0000 
7.0000 
7.0000 
7. 0000 
7.0000 
7 .0003 
7.0000 
7.0000 
7.0000 
7.0000 
7.0000 
7.0000 
7.0000 
7 .0003 
7.0000 
7.0000 
7 .0003 
7.0000 


1 .0000 

C.  99 ?4 

0 . o « 8 n 

0.9  8 56 

0 . 9 8 2 3 

0.9789 

C.  9748 

0.9692 

0.9590 

0.  9097 

0.7316 

0.4  330 

0.4300 

0.4300 

0.4  300 

0.4300 

0.4380 

0.4300 

0.4300 

0.4  330 

3.430  0 

0.4300 

0.  430  0 

0.4300 

0.4300 

0.4  3 00 

0.4300 

0.4300 

0.  430  0 

0.4300 

0.4300 


0 .33 OP-1  3 
9. 4118-10 
0. 45 5 E- 10 
3 . 499F-1  0 
3.  547F-10 
0 . 60 4 E- 10 
0. 679E-10 
0.  796E-10 
0 . 1.06F-39 
9.  4 ^ 9 F - 0 9 
0. 781 E-07 
C. 674E-03 
0.  674-=- 93 
0 . 674E-03 
0. 674c-03 
. 6748-03 
, 6T4F-T13 
674E-03 
0. 674E-03 
0 . 674E-03 
O'. '67477- 03“ 
0. 674  8-0  3 
0.  674^-03 
0.  674E-J3 
0 . 6748-33 
0. 674E-03 
0. 674E-03 
0. 674E-03 
0.  674!'- 03 
0.6748-03 
0. 674E-03 


0, 

0 

n. 


3 .0 

0. 1 91  F-07 
0 .407 E-07 
3 .6+68-0 7 
0.  907? -07 
3 .1196-06 
0. 151 F-06 
0. 188F-06 
0 .233  F-0  6" 
0.331 F-06 
0 .5018-05 
3.1878  00 

n. 00 

0.860F  00 
0 . 123F  01 
0.  153E  01 
0.1378  91 
0.22  IF  01 
0 .2548  01 
0 .2338  01 
'0Y7Z2E  01“ 
0.355E  01 
3.339F  01 
0.4238  01 
0 . 4 5 6 r 01 
0.4  906  01 
0.524F  01 
0.557E  01 

0. 59rE“BI 

9 . 6 2 5C  01 
3.6599  01 


0.0 ' 

0.617E-39 
0.494E-08 
0 .167E-07 
"'TTVB'OTF-OT"" 
0.  7728-07 
0.133  E— 0 6 
0.212E-06 
■■~tr.'3T6F=inr- 
0.4  5 OF- 06 
0.61  78-06 
0.8  2 IE- 3 6 
"~OTT06F-CrT 
0.  13  5E-0  5 
0. 169F-05 
0.208E-05 
0.2  5 28-0  5 
0.301 F-05 
0.357E-05 
0 .41 9F-0  5 

0.563E-05 
0.646E-05 
0.  735E-05 

0.93  6 F-0  5 
0.105E-04 
0 . 1 1 7 r-  0 4 
0.72  98-04 
0.  143E- 04 
0.157E-04 
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In  Table  7 the  probabilities  of  first  aircraft  failure 
for  a single  airplane  and  in  a fleet  of  50  airplanes  are  shown 
for  the  case  when  no  inspections  are  performed. 

TABLE  7.  PROBABILITY  OF  FIRST  AIRCRAFT  FAILURE  (NO  INSPECTIONS) 


i,T,ptn » tp{  i» 

Single 

. j__^ 

Fleet 

l 

500.0 

0 . 15  522E- 07 

0.826  OFF- 05 

2 

1C00.  0 

0.33  043E-07 

0. 1652 1F-05 

3 

1 50C.  0 

0.  49586E-  07 

0. 247  93E-05 

4 

2000.0 

0.66197E-07 

0.  33098F-05 

5 

2500.0 

TT/srowE-or 

0.41 49 5 E -05 

6 

3 000.0 

0.99535E-  07 

0.  49767F-05  ' 

7 

3500.0 

0 . 11 640 E-  06 

0 . 58.2 0 QE-05 

8 

4000.  0 

0.1345  8E-06 

0 . 672  8 8E-05 

9 

4500. 0 

0.155  30E  - 06 

0.TT649E-TI5F 

10 

5000.0 

0 .17 369E-06 

0 . 8684  5E-05 

11 

5500.  0 

0.  1941  8E-  06 

0.  9709  OF -05 

12 

6000.0 

0.25  7 74E-  06 

0.12887F-04 

13 

6500.0 

0 . 73939F-  06 

0 . 365T9F-04 

14 

7000.  0 

0.  27  1 08c-  05 

0.  13553E-03 

15 

7500.  0 

0.83435E-05 

0.41709E-03 

1 6 

8COO.  0 

0.21  126F-04 

0 . 1055  8E-02 

17 

~8500.0 

0 .45  866r-04- 

0.  2290  7F-02 

18 

9000.0 

0.39012E-04 

0.44409E-02 

19 

9500.  0 

0.  15871E-03 

0.79041E-02 

20  ] 

.0000.0 

0 .26  46  IE  - 03 

0.  13145E-01 

21 

0500.0 

0 .4T823E-U3 

13 .2(159  9~E-DT 

22  1 

.1000.  0 

0.63  304E-  03 

0.  3116  6E-01 

23  ] 

L 1 500.  0 

0.92  38 1 E-  03 

0.4516  IE-  01 

24  1 

2000.  0 

0.13  0 75E-  02 

0.63325E-01 

25  ! 

,2500.0 

0.18(F27F-'U2~ 

0 • 862  7 CE  — 01 

25 

3000.0 

0 .24298E-02 

0. 11453E  00 

27 

L 3 500,  0 

0.  32106E-  02 

0.14853E  00 

28 

,4000.0 

0.41690E- 02 

0 . 188  5 IE  00 

29 

.4500.  0 

CT.  53304F-U2 

0 *234-5  GE  00  i 

30 

15000.  0 

0.  672 16E- 02 

0. 2862 4E  00 

The  probabilities  of  first  aircraft  failure  for  a single  airplane 
and  a fleet  of  50  airplanes  with  5 scheduled  inspections  are 
given  in  Table  8 . 
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TABLE  8 


PROBABILITY  OF  FIRST  AIRCRAFT  FAILURE 
(INSPECTIONS  INCLUDED) 


I , T,  P ( I 


1 

2 

3 

4 

5 

6 


* » TP ( I ) Single  Fleet 

2500.  C 0 . 8?  990E-  07  0.41495E-05 

5000.0  0.17184=-  06  0.  8591  9E-05 

750 0.  0 0 .72 494E-06  0. 36246E-04 

10000.  0 0.  10433E-04  0.521  50E-03 

12500.0  0.53212E-94  0.  26571E-02 

15000.  0 0 .1 6 8 89 E -03 0.84099E-02 
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SECTION  4 


AIRCRAFT  RELIABILITY  WHEN  INITIAL  CRACK 
SIZE  IS  A RANDOM  VARIABLE  (AREL2 ) 

The  purpose  of  AREL2  program  is  to  determine  the  proba- 
bility of  first  aircraft  failure  either  of  a single  airplane  or 
of  an  airplane  in  a fleet.  The  probability  of  aircraft  failure 
is  constructed  on  the  condition  that  the  initial  fatigue  crack 
size  is  a random  variable  for  which  the  probability  density 
function  is  obtained  by  appropriate  transformation  of  the  Wei- 
bull  distribution  of  the  time  to  fatigue  crack  initiation.  The 
probability  of  failure  model  allows  to  include  rigorous  base 
type  inspections  which  are  performed  at  equal  time  intervals 
over  the  design  life  of  the  aircraft  and  cursory  inspections 
which  are  performed  at  specified  time  intervals  between  the 
rigorous  inspections.  The  probabilities  of  aircraft  failure 
are  calculated  at  rigorous  inspections.  The  program  includes 
options  to  study  the  effect  of  proof  test  and  input  load  level 
truncation  on  the  probability  of  first  aircraft  failure.  A 
simplified  flow  chart  of  the  program  is  shown  in  Figure  3.  The 
following  paragraphs  discuss  the  input  data,  the  structure  of 
the  program,  and  the  output. 
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4 . 1 INPUT  DATA 


There  is  one  input  data  card  which  is  of  the  same  format 
as  shown  in  Table  6 except  the  parameter  XY  now  ranges  from  1 
to  34.  The  remaining  data  parameters  are  specified  at  the  be- 
ginning of  the  main  program.  The  main  contents  of  the  data 
parameters  for  AREL2  program  is  the  same  as  for  ARELl  program. 
These  are  described  in  Sec.  3.1.  The  additional  data  parameters 
are  as  follow:  NBN  = number  of  subdivisions  on  the  probability 

density  function  G(A0)  corresponding  to  initial  fatigue  cracks; 
AL  = lower  limit  of  AO  for  G(A0);  AU  = upper  limit  of  AO  for 
G(A0);  ACT  = truncation  level  due  to  proof  test  on  G(A0),  CSRY  = 
0.0  implies  no  cursory  inspections,  CSRY  = 1.0  implies  cursory 
inspections  included;  NR  = number  of  rigorous  inspections;  NC  = 
number  of  cursory  inspections;  CA1  = minimum  crack  size  which  is 
detected  by  cursory  inspection. 


4 . 2 PROGRAM  DESCRIPTION 

The  Aircraft  Reliability  Program  AREL2  consists  of  a 
main  program  unit  and  eleven  subroutines  DENST,  HCOM,  HSIN,  QSF, 
CIN,  HIN,  ERF,  OMEXP,  FF  and  FC  which  are  described  in  the 
following  paragraphs.  The  program  listing  is  given  in  Appendix  C. 

4.2.1  The  Main  Program  Unit  AREL2 
The  main  program  reads  and  initializes  data  as  described 
in  Sec.  4.1.  The  first  part  of  the  program  is  designated  for 
the  sensitivity  study  of  the  34  parameters  in  the  reliability 
model.  This  is  accomplished  by  specifying  the  value  of  XY  from 
1 to  34.  Each  number  corresponds  to  a particular  parameter  in 
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the  system  which  are  described  in  Sec.  4.1.  The  combined 
input  load  Q is  determined  in  the  same  fashion  as  described 
in  Sec.  3.2.1.  Then  with  a call  to  subroutine  DENST  a 
probability  density  function  G(A0)  for  the  initial  crack 
size  AO  is  generated.  The  crack  size  A(I) , the  residual 
strength  R(I)  and  the  failure  rate  H(I)  are  calculated  in 
terms  of  AO  in  a similar  fashion  as  described  in  Sec.  3.2.1. 
The  summation  of  the  failure  rate  H(I)  is  obtained  by  call- 
ing subroutine  QSF.  The  probability  of  first  aircraft  fail- 
ure with  rigorous  and  cursory  inspections  included  is  cal- 
culated from  subroutine  CIN.  Similarly,  the  probability  of 
first  aircraft  failure  where  no  inspections  or  only  rigor- 
ous inspections  are  performed  is  determined  by  calling 
subroutine  HIN.  The  effect  of  proof  test  on  the  probabil- 
ity of  failure  is  determined  by  specifying  ACT  to  be  equal 
or  less  than  AU. 

4.2,2  Subroutine  DENST  (G,  NN,  AU,  AL,  DA) 

Subroutine  DENST  calculates  the  probability  density 
function  G(A0)  of  the  initial  crack  size  AO.  This  is  achieved 
by  transforming  the  probability  density  function  W(T)  corre- 
sponding to  the  time  to  crack  initiation.  The  output  of  the 
transformation  is  stored  in  G.  The  other  arguments  of  the 
subroutine  are 

NN  = number  of  division  points  on  W and  number  of 
output  points  on  G; 

AU  = maximum  specified  initial  crack  size  value; 

AL  = lower  limit  crack  size  of  G. 


-32- 


4.2.3  Subroutine  HCOM  (H,  ARO , AR1 , . . 


•r 


AR8 , P,  AN) 


Same  as  described  in  Sec.  3.2.3. 

4.2.4  Subroutine  HSIN  (H,  ARO,  ARl , ■ . ♦,  AR8 , P,  AN,  UPC) 

Same  as  described  in  Sec.  3.2.4. 

4.2.5  Subroutine  QSF  (H,  Y,  Z,  NDIM) 

Same  as  described  in  Sec.  3.2.5. 

4.2.6  Subroutine  CIN  (A,  H,  G,  CA1,  NC,  NR,  NN,  ITO) 

Subroutine  CIN  determines  probabilities  of  first  aircraft 
failure  with  both  rigorous  and  cursory  inspections  included.  The 
arguments  of  the  subroutine  are 

A = input  data  - fatigue  crack  size  which  is  calculated 
in  the  main  program; 

H = summation  of  the  failure  rate; 

G = probability  density  function  of  the  initial  crack 
size  which  is  determined  by  subroutine  DENST; 

CA1  = minimum  crack  size  which  is  detected  by  cursory 
inspection  and  is  specified  as  input; 

NC  = number  of  cursory  inspections  in  an  interval  between 
rigorous  inspections; 

NR  = number  of  rigorous  inspections; 

NN  = total  number  of  inspections  (cursory  and  rigorous) ; 
ITO  = an  index  for  output  print. 

4.2.7  Subroutine  HIN  (A,  H,  G,  NN,  ITO) 

Subroutine  HIN  calculates  probabilities  of  first  aircraft 
failure  when  no  inspections  or  only  rigorous  inspections  are 
performed.  The  meaning  of  the  arguments  for  this  subroutine 
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are  the  same  as  for  subroutine  CIN. 


4.2.8  Function  ERF (X) 

Function  ERF(X)  is  the  same  as  in  AREL1  program. 

4.2.9  Function  OMEXP(X) 

Function  OMEXP(X)  evaluates  OMEXP(X)  = 1 - EXP ( X) 
utilizing  an  expansion  procedure. 

4.2.10  Function  FF(X) 

Function  FF(X)  calculates  the  crack  detection  probability 
corresponding  to  rigorous  type  inspections.  The  necessary  data 
for  this  function  is  specified  in  the  main  program. 

4.2.11  Function  FC(X) 

Function  FC(X)  specifies  the  crack  detection  probability 
for  cursory  inspection.  The  input  argument  in  terms  of  the 
crack  detection  size  is  specified  in  the  main  program. 

4 . 3 PROGRAM  OUTPUT 

The  AREL2  program  generates  printed  output.  The  input 
data  XY,  YZ1,  YZ2,  YZ3  is  printed  for  a chosen  example.  The 
number  of  rigorous  inspections  NN,  the  crack  size  A,  the  failure 
rate  H,  the  probability  of  first  aircraft  failure  of  a single 
airplane  PA,  and  the  probability  of  first  aircraft  failure 
TPA  in  a fleet  of  airplanes  are  printed  as  shown  in  Tables  9 
and  10  corresponding  to  fail-safe  and  slow  crack  propagation 
models,  respectively. 
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TABLE  9.  CRACK  SIZE  A,  FAILURE  RATE  H,  PROBABILITY  OF  FAILURE 
FOR  SINGLE  AIRPLANE  PA,  PROBABILITY  OF  FAILURE  IN  A 
FLEET  TPA  (FAIL-SAFE) 


HN='"L4 

r a 

11  0 . 42  704  66E-0 1 
21  0.46673  49E-0 1 


3l  0.53641  76E-01 
41  0 . 75  46765E-0 1 
51  0. 15  828  03  E 01 
61  C»  7000000E  01 


0. 4072061  E-0  7 
0.9069237  E-0  7 


O.700O0OOE  01 
0.7000000  E 01 
0.7000900E  01 
000000  E 01 


0.16  i’3  89b  E-06 
0.23  3108  9E- 06 
0 .5295597  E-05 
0.52 07602  E OC 


C.  3394086E 
0.68155  83F 


PA 

-07  0. 
•07  0. 


0. 1024292b -06  07 
0.  13671  63E-06  0. 
0. 171G062E-06  0. 
0.  2 053  084 E -06  0. 


TPA 

1697041 E-05 
3407785E -05 
51? 1445 E-05 
6 83  5792  E-05 
8 55  0273E - 05 
1 026537  E-04 


71 

81 

91 

0\ 


.11 

121 

131 

JL41 


3$ 


151 


O'OOOCOE  01 
0.7000000E  01 
0.  7000000  E 01 
0. 7000000F  0 1 


0 .1154295  fc  01 
0.1867 82 7E  01 
0.2541360F  01 
0 .32  1 489'2-E  01 


0 . 70000 00 E 01 


U7MW5TTT 
0.4561957E  01 
0 .52  35490  E 01 
0.  5909C22  E 01 
0 .65  82555  E 01 


U7Z39'6b  73F-06  0. 
0.2  7404  77  E-06  0. 
0. 3 084521 E-06  0. 
0.34287  97E-06  0. 


WO. 
•06  0. 
-06  0. 
-06  0. 


1198?  792-04” 
1 37  022  9E -04 
1542243  E-04 
1 71438  3E - 04 


T). '37  731  34E 
0.  4117510E 
0.4461873E 
0.48062  33 E 
7T75T5U 5 94F-06  0.  257 526 3t - 04 


T"8  8 5 5 5 0 E - 0 4 
2C5  8732E-04 
223091  OE-04 
2403087E-04 


TABLE  10.  CRACK  SIZE  A,  FAILURE  RATE  H,  PROBABILITY  OF 

FAILURE  FOR  SINGLE  AIRPLANE  PA,  PROBABILITY  OF 
FAILURE  IN  A FLEET  TPA  (SLOW  CRACK  GROWTH) 


<N«  14 

T , A H PA  T PA 

11  0.4270466 E-01  0 .3304082 E-0 7 0 . 33040 45E-0 7 0.  1652020E-05 

21 DT46 673  49E-01  0.6608151  E-07  0.660/9  74^-07  0.  3303 98TE=TT5” 

31  0.5364176E-0  1 0 .99  1 2 242 F- 0 7 0.  9 9 1 19  18E-0 7 0.  4 95 5 946E - 05 
41  0.7546765 E-0  1 0 .132 1632 E-06  0 . 1 32 1 5 84E-06  0. 6607897E- 05 
51  0.15  828  03E  01  0.  I 76591  0E-04  0.  1 651  9 74  E-06  0.  825 983  8 E-05 

61  O.IOOOOOOE  04  0.7670562F  03  0.  19840  70E-06  0.  9920299E- 05 

— 7T  0.  tOO 0 0 COE  04  0 . 1 76rf-ffi5  E 04  0 . ? 322 5 12  E-06  0 . 1 16 1 249E -04 

81  O.IOOOOOOE  04  0 .2767054E  04  0.  2671  l 56F-06  0.  1 33  5569E-04 

91  O.IOOOOOOE  04  0.3767052  E 04  0 . 30  2 1 1 81 E - 06  0.  1 5 1 0579E- 04 

101  0. 1 0000  COE  04  0.4767 03 9E  04  0 . 3 37? 5 25 E-06  0 . 1686247 E-04 
III  O.IOOOOOOE  04  0.5767  02  OE  04  0.  3 7?  5 5 74E- 06  0.  1 862  76  9E-04 

121  O.IOOOOOOE  04  0 .6767000E  04  0 . 4078  3 62  E-06  0 . 203  91 60E  - 04 

131  O.IOOOOOOE  04  0.7766980E  04  0. 4431 927 E-06  0. 221 5933 E-04 

141'  O.IOOOOOOE  04  0.8766961E  04  0 .478  5 2 539^05  0.-239  26 172  - 04 

151  O.IOOOOOOE  04  0.9766941 E 04  0 . 5 1 388 96E-06  0 . 2 5694 1 4E-04 
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APPENDIX  A 


MANUVER  LOADING  . 

DIMENSION  A (3276 8),  INV(4096)  , S ( 4096 ) , BETA ( 6 l 
DIMENSION  A1(8),A2(8),A3(8J,A35(8),A4(8),A4  5(8),A5(8) 
DIMENSION  Bl (8) » E2( 8 ) , B3( 8 ) , B35C  8) , B4 ( 8) ,84  5( 8) , 95 ( 8 ) 

01  MENS  ION  XI (6) , X2(6),X3(6) , X3 5 ( 6) , X4 ( 6 ) ,X45(6) ,X5(6) ,XX(6,5) 
DIMENSION  Yl { 6 ) , Y2 ( 6 ) , Y3 ( 6 ) # Y35 ( 6) , Y4 ( 6 ) , Y45 ( 6J * Y5 { 6 ) , YY ( 6,5 } 
DIMENSION  NPK(8I ,MPK(8 ) , NPEAK( 6) ,MPEAK(6) 

DIMENSION  AA(8, 5 > , BB (8 , 5 ) , AN ( 5 I , BN ( 5 ) 

COMMON  /AXC/ZI200)  

CATA  BETA/O. 0,0. 25, 0.5  ,0.75,0.8,0.9/ 

NSAMPL  — NO  OF  SAMPLES  USED  FOR  EACH  8ET A VALUE 

NBETA  — MO  OF  BETA  VALUES 

WUP  — UPPER  BOUND  FREQUENCY  IN  HERTZ 

WL  — LOWER  POUND  FREQUENCY  IN  HERTZ 

N — NO  CF  DIVISIONS  IN  SPECTRAL  DENSITY 

NPT  — NO  DF  SIMULATED  TIME  POINTS 

Ml  --  2**M1=NPT 

NSIGN  — 0 USE  WHITE  NOISE  SPECTRUM 

1 USE  USER  SUPPLY  SPECTRUM 
TIME  --  TOTAL  SIMULATED  PERIOD  IN  SECONDS 
REOUCE  — REDUCTION  FACTOR  OF  THE  NEGATIVE  PART  OF  THE 
TIME  HISTORY 

REDUCE=0. 25 
PI=3. 141593 
NSAMPL=  8 
NBET  A=1 
WUP=  0. 2 
N=  2048 
NPT= 16384 
Ml=14 
NS  IGN= 1 

OT=FLOAT (N ) / ( WU  P *EL  0 AT (NPT I ) 

WL=0. 007 
TIME=DT*NPT 

C CALL  AREA  TO  FIND  THE  AREA  OF  SPECTRAL  DENSITY  IN  ORDER  TO 

C NORMALIZE  S(W> 

CALL  AREA  (S,512  ,WUP,WL> 

A(  11=100000. 

A ( 2 ) =0. 

A ( 3 ) * • 2 
SCAL E=1.0 

C CALL  GRAPH  (A, S, 512, 2,'  ',2,*  ',2,'  »,2,«  •, SCALE) 

WR IT E( 6 , 31 2 ) TIME, DT, WUP, WL 

312  FORM AT ( l OX ,' PER  1 0D= ' ,F10.2,4X, *DT=» ,F6.4,4X, *WU  = • ,F5.3, 

1 4X, * WL= ' ,F5. 3/  ) 

NP=4Q0 

C NPOTS  — NO  OF  POINTS  IN  EACH  GRAPH 

C NPLOTS  — NO  OF  GRAPHS  TO  BE  PLOTTED 

NPOT S=1025 
IX=3234567 

DO  100  1=1, NBETA  _ 

WU=WUP 
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00  80  J = i » NSAMPL . ... _ 

CALL  SIMULT  ! A, N ,NPT ,WU, WL , IX, Ml ,NS IGN, INV , S) 

00  72  KK=1,NPT 

IFIA!KK).LT.0,01__A!KK)=A!KK)*REDUCE  __  _ 

72  CONTINUE 
WR  ITE!  6 » 20  ) 

20  F0RMAUJLH1)  __  __  _ 

IF  ! J.GT.l)  GO  TO  15 

C WRITE  OUT  PART  OF  THE  SAMPLE  TIME  HISTORY  (1ST  SAMPLE) 

NS=NPT/8 _ _ __  __  

WRITE!  6,  74  ) (All  I), 11  = 1, NS) 

74  FORMAT! ( IX, 24F5. 2 ) ) 

CALL  PLOT  !NP,0T,A,0.0,1) 

15  CALL  PICK  !A  »NPT,S»NPKC J) ) 

CALL  STDEVI(1,NPK!  J ) , S ,STD,  RMS , AH  J),A21  J)  ,A31  J),A35iJ) 

1 ,A41J)  ,A45(  J),A5!J),AN) _ 

DO  34  L=1 , 5 
34  AA! J,L)=AN!L) 

WRITE! 6, 25)  BETA ! I) , NPK! J ) , WU, WL  ,STO,RMS 
25  FORMAT ( //5 X, • BETA=*  »F6.3,5X, *NPEAK= • , 16, 5X, *WU=* 

1  F8.4,5X,*WL=*,F8.4,5X,*  STOE • , E14.7 , 5X , • RM S= • , E14. 7/ ) 

WRITE !6, 30)  Al!J),A2(J)«A3!J)»A35!J),A41J)tA45lJ),A5!J), 

1 I AA ! J, L ) ,L=1 » 5 ) 

30  FORMAT (5X,* MEAN  1 = *, F7.3, 3X, *MEAN2=* ,F7.3, 3 X, 'MEAN3  = ' ,F7.3, 
1 3X,  «MEAN35=»  ,F7.3,3 X,*M£AN4=*  ,F7. 3,3 X,  ' MEA N45=J  , F 8_*  3 _ 

2 3X»  * MEAN5  = * , F8 . 3/5X , • MEAN6= », E 1 2. 5 , 2X, • MEAN7=* , E12. 5,2X, 

3 *MEAN8=*,E12.5,2X,*MEAN9=* ,E12.5,2X, *MEAN10=*,E12.5/) 

WRITE  !6, 70  ) TAIL ),L=i,NP) . _ 

70  F0RMAT!/10(2X,F10*5) ) 

WR  IT  E 1 6 , 70 ) !S(L),L  = 1,50) 

_ CALL  PICKM_(  A,NPTji^,JMPK!  J ) ) 

CALL  STDEVI ! 1,MPK1  J)  , S ,ST0,RMS  , Bl!  J ) , B2  ( J)  , 63 1 J ) , B35  ( J ) , 

1 B4!  J ) , B45! J ) , B5 ! J ) , BN  ) 

00_  44  L=I , 5 _ _ _ __  . _ _ 

44  B B I J , L ) = BN! L ) 

WRITE! 6,25)  BETA ! I ) , MPK( J ) , WU, WL  , STD, RMS 

W R IT E ! 6 , 30  ) „B1 ! J ) , B2 ( J T,  B3!  J ),  B35!  Ji,  B4i  J ) , B45!  J ) , B5  ! J i _ 
l , ! BB ( J , L ) , L = 1 , 5 ) 

C WR IT  E! 6, 70  ) ! S! L ) ,L=1, 50 ) 

CALL  STATJS  I A,NPT,STDV,RMSV)  _ 

C TO  FIND  THE  STANDARD  DEVIATION  OF  THE  PROCESS 
WRITE!6»212 ) STDV,RMSV 

212  FORMAT ! //5X , 1 STANDARD  DEVIATION  OF  THE  TIME  HISTORY* 

1 //5 X, * STDEV=* , E14.7 ,5X»*R.M.S.='» E14. 7 ) 

80  CONTINUE  . 

NSAM=NSAMPL  _ _ „ 

CALL  MEANT  ( N SAM , A1 , A2 , A3 ,A 3 5 , A4 , A45, 45  , AMI , AM2 , AM3 , 
l AM35»AM4»  AM45, AM5,AA, AN ) 

JrfRITE<6,95_)  AMI  , A M2  , AM3,  AM3  5 ,AM4, AM45,  AM5,  1 ANIL  ) ,L=1,5)  _ 
9 5 FORMAT 1////10X, * AM  • , 6 1 3X , E 14. 7 ) / 1 5 X , 6 ! 3 X, E 1 4. 7)7) 

XI!  I )=AM1 

X2J I )=AM2  _ ______  _ _ 

X3  ( I )=AM3 
X35!  I )=A  M35 

X4!I)=AM4  ______ 

X45! I )=AM45 
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X 5 ( I ) = AM5  

DO  82  L=1 » 5 
82  XXI I,L)=AN!L) 

CALL  MEAN7  1NSAM  ,B1 , B2  , B3,  B35,  B4,  B45,  B 5,  BM1 , BM2 , BM3,  __ 
1 BM35, BM4, BM45, BM5, BB»  BN ) 

WR  IT  E ( 6,  96  ) BM1 , BM2,  BM  3,  BM35 ,BM4,BM45,BM5,  ( BNIL) ,L=1,5) 

96  FORMAT (//l OX,' BM  » ,6< 3X , E 14.7 ) / 1 5X ,6 1 3X , E14. 7 ) / ) 

YiU)  = BMl 
Y2C l)=BM2 

Y3(  I ) = BM3  . . _ _ 

Y35( I ) = BM35 
Y4( I ) = BM4 

Y451I)=BM45  _ _ ______  _ 

Y5 ( I ) = BM5 
DO  81  L = l,5 

81  YYI l,L»=BN!L)  _ _ _ 

NPEAK (I ) =0 
MPEAKl I I =0 

DO  98  I I = 1_*NSAM  ___  __ 

NPEAK! I )=NPEAK( I )+NPK( II ) 

98  MPEAKU  )=MPEAK(I)+MPK(  II  I 

NPEAKII  )=NPEAK(  I )/FLOAT(NSAMJ _ 

MPEAKU  )>MPEAK(  I J/FLOATINSAM  ) 

100  CONTINUE 

WR  ITE(  6, 105 1 _ __  _ 

105  FORMAT (1H1////42X, 'CASE  I'// I 
WRITE! 6, 106) 

106  F0RMATC30X,  • AVERAGE  OF  POWER  OF  RI$E_AND  FALL*/  _ 

1 /30X » • AND  NUMBER  OF  RISE  AND  FALL  * //43X , * POWER’ / ) 

WRITE(6, 110 ) 

110  FORMAT  !15X,'1.0'  ,5X,  '2.0', 5X,' 3.  O'  , 5X»  '3.  5 ' j 5X , *4.  0 _ 

1 5X, '4.5' , 5X, *5. 0* , 5 X, 'RISE  £ FALL • /6X, • BETA • ) 

DO  120  1=1 , NBETA 

120  WRITE  !6 , 140 ) BETA!  1 1 1 X 1(  I J , X2(  IJ„,  X3  1 1)  , X35I  I JjlX4U  ) t 

1 X45 II), X5 ( I ) , NPEAK ( I ) 

140  F0RMAT(/6X,F4.2,2X,7F8.3,2X, 17 > 

141  F0RMAT!/6X,F4.2»2X,51F10.3,lXi,2X,  16)  _ 

WRITE! 6, 145) 

145  F0RMAT!1H1////42X,'CASE  II'//) 

WR  ITE!  6,106)  __  - _ _ 

WRITE  16,110) 

DO  150  1=1,  NBETA 

150  WRITE!  6, 140 ) BE  T A1J  » , Y 1 1 I ) , Y2(  I ) , Y3  1 Li*  Y35  ( I ) , Y4  1 1 )_  

1 , Y45!  I ) ,Y5!  I), MPEAKU  ) 

WR  ITE16, 105 ) 

WR  ITE!  6, 106 ) _ _ 

WRITE  16, 111 ) 

111  FORMAT  1 1 5X, '6.0'  , 8X,  '7.0 ' , 8 X , • 8. 0 • , 8X , » 9. O' , 8X,  ' 10. • , 

_1  8X» 'PEAKS'/ 6X, 'BETA')  _ 

DO  121  1=1, NBETA 

121  WR I TE 1 6, 141 ) BETA1I),!XX!I,L),L=1, 5), NPEAKII) 

WRITE!6,145)  __  _ _ _ 

WRITE (6, 106) 

WRITE(6,111 ) 

DO  151  1=1, NBETA  _ _ _ 

151  WRITE(6,141)  BETA  1 1 ) , ( YY 1 1 , L ) ,L=1 , 5) ,MPEAK 1 1 ) 
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STDEV=1.  - 

CALL  PICKMX  (A*NPT*S(1)»LLL*S(2049)»MMM) 
WRITE 16, 123)  LLL  ,ISI L>  ,L=1,LLL  ) 


WRITE! 6tl24 ) MMM , ISll+2048) , L=1*MMM) 

123  FORMAT  ( 1H1  /10X,  • NO.  OF  PEAKS=’  , I 6/  ( 2X , 10F8.  4 )) 

124  FORMAT ( IH1 /10X, ' NO.  OF  TR0UGHS=’ , I6/f 2X, 10F8.4) ) 

RATIO  — 1 GIVE  TOTAL  CROSSINGS  OF  THE  HISTORY 

~ 3600/ PERIOD  GIVE  CROSSINGS  PER  HOUR 


RATI0=1. 
CALL  CROSS 


(STDEV»S II) i$C  2045!  ♦LLL ,NMM, VZERO, RAT  10 ) 


RAT I0=3600« /TIME  . 

CALL  CROSS  I ST DEV, Sill ,S I 2049) *LLL  ,MMM»  VZERO  t RAT  10 ) 

S_TQP  — - 

END 


SUBROUTINE  STDEVIINST.NPT.^UN, STD,RMS,A1,A2,A3, A35,A4,A45,A5,AN) 

DIMENSION  FUNC1 ) ,AN! D 

ASUM1=0.0 

ASUM2=0.0 

ASUM3=Q.0_ 

ASUM35=0.0 
ASUM4=0.0 
ASUM45=0.0 
ASUM5=0.0 
DO  10  1=1,5 
10  AN( I )=0. 0 

XPO INTaFLOATCNPT-NST  *1 ) 

DO  5 JT=NST ,NPT 

AA2=FUN I JT I *FUN I JT ! _ _ 

AA3=FUN( JT)*AA2 
AA35=FUNIJT>**3.5 
A A4=  AA2*AA2 
AA45=FUN(JT)  **4.5 
A A5=AA3*AA2 

ASUMDASUM1+FUNI  JT)  _ 

ASUM2=ASUM2+AA2 
ASUM3=ASUM3*AA3 
A SUM 35  = ASUM3  5+A  A 35 
ASUM4=ASUM4*AA4 
ASUM45=ASUM45+AA45 
ASUM5=ASUM5+AA5 
ANI1 )=AN II )+AA3*AA3 
AN! 2 ) = AN 1 2 ) +AA3  5*AA  35 
ANI3)=AN(3)+AA4*AA4 
ANI4)=AN(4)*AA45*AA45 
AN! 5 )= ANI 5 ) + AA5*AA5 
5 CONTINUE 

ADA  SUM1/XP0INT 

A2=ASUM2/XP0 I NT 

A3=ASUM3/XP0INT 

A35-ASUM35/XP0INT 

A4=ASUM4/XP0 INT 

A45=ASUM45/XP0INT 

A5=ASUM5/XP0INT 
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00  15  K=l » 5 
15  AN(K)=AN(K) /X POINT 
BMS=SQRT ( A2 ) 
STD=SQRT(A2-A1*A1) 
RETURN 
END 


SUB ROUT INE  MEANT  ( NSAMPL , A 1 , A 2 , A3 , A35 , A 4 , A45 , A5 , AMI  TAM  2 , 
l AM3,AM35,AM4,AM45,AM5,A.A,AN) 

Of  MENS  I ON  Al(l)  , A2(l ), A3 ( 1 ) , A35 { 1 ) , A4 ( 1 ) ,A45(1) ,A5(1 ) 
DIMENSION  A A ( 8 , 5 ) , AN  ( 5 ) 

AM1=0. 

AM2=0. 

AM  3=0. 

AM35=0. 

AM4=0. 

A M45=0. 

AM5=0. 

DO  10  1=1,5 
10  AN ( I )=0. 

DO  90  K=l, NSAMPL 
AM1=AM1+A1 (K ) 

AM2=AM2+A2(K> 

AM3=AM3+A3 I K ) 

AM 35= A M3  5+A35 ( K ) 

AM4=AM4+A4IK ) 

AM45=AM45+A45IK ) _ 

AM5= AM5+A5 ( K ) 

DO  80  L = 1 » 5 

80  AN(L)  = AN(LH-AA(K,L) 

90  CONTINUE 

S A MPL= FLOAT! NSAMPL ) 

AM1=AM1/SAMPL 
AM2=AM2/SAMPL 
AM3=AM3/ SAMPL 
AM35=AM3  5/SAMPL 
AM4=AM4/SAMPL 
AM45=AM45/SAMPL 
AM5=AM5/SAMPL 
DO  85  L=1 , 5 
85  AN(L)=AN(L)/SAMPL 
RETURN 
END 
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SUBROUTINE  PICK  ( fit NP, PEAK, NPK ) 


PICK  UP  DIFFERNECE  BETWEEN  SUCESSIVE  RISE  AND  FALL 
A — INPUT  ARRAY  FROM  WHICH  RISE  AND  FALL  ARE  TO  BE  FOUND 

NP — NO  OF  PQJNTS  IN  LNPUTJARRAY  A __ 

PEAK  — ARRAY  TO  STORE  THE  RESULT  OF  i PEAK ( I ) -TROUGH 1 1-1  I ) 
NPK  — NUMBER  OF  VALUE  FOUND 


DIMENSION  A(  1 ) t PEAK II) 

PEAK  (1)  = A(l) 

I FLAG*  l 

IF(A(1)»GT*A(2)  ) IFLAG=0 
I PK=2 

D0_100  I = i ♦ NP _ 

IFUFLAG.LT.1)  GO  TO  6 0 
PEAK  TESTING 

IFLA<  n.LT,A(JdJ)_GO„TO  _100 __ 

PEAK  I IPK )=A I II 

PEAK(IPK-U  = ABS(PEAKC  IPK  )-PEAK  { IPK-1)  ) 

__  I PK*  I PK+  I _ 

IFLAG=0 
GO  TO  100 

TROUGH  TESTING 

60  IFTA(ii.GT.AUU))  GO  TO  100 
PEAK ( IPK I =A  f I ) 

PEAK!  IPK-1  ) = ABSJPEAK(JPK)-PEAKJJPKrlJJ 

IPK* IPK+1 

IFLAG=1 

100  CONTINUE 

NPK* IPK-2 
RETURN 

END  _ _ _ _ 


SUBROUTINE  PICKM  I A , NP  .PEAK, MPK ) 

PICK  UP  THE  DIFFERENCES  BETWEEN  SUCESSIVE  MAXIMUM  RISE  AND 
MINIMUM  FALL  < 2ND  METHOD  OF  PICKING) 

A " INPUT  ARRAY  FROM  WHICH  RISE  AND  FALL  ARE  TO  BE  PICKED 

NP  — NUMBER  OF  POINTS  OF  INPUT  ARRAY  A 

PEAK  _r-  ARRAY  USE  _T0 _H.OLD.THE  RISE  AND  FALJL  DIFFERENCES 

MPK  — NUMBER  OF  RISE  AND  FALL  DIFFERNENCE  FOUND 

DIMENSION  Ajll.PEAK (1 ) _ ______ 

AMAX=0. 

AMIN=0. 

LFLAG_=1  _ ___  __ 

IFCACll.LT.O.)  I FLAG =0 
I PK=1 
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00  100  1=1, NP  _ 

IFUFLAG.EQ.O)  GO  TO  40 
IF(A(  D.GT.AMAX)  AMAX=A(I) 
I F ( A I I +1 ) .LT #0#  ) GO  TO  50 
GO  TO  100 

50  IFL4G=0 

PEAK(IPK)=ABS<AMAX-AMIN) 

I PK= I PK  + 1 
AMIN=1.0 
GO  TO  100 

40  IF(A(  I)  .LT.AMIM  j AMIN=A(I> 
I F f A C I + l ) • GT .0. ) GO  TO  70 
GO  TO  100 
70  I FLAG=1 

PEAK(IPK)=ABS( AMAX-AMIN ) 

1 PK= I PK+ 1 
A MAX =-1.0 

100  CONTINUE 

MPK= I PK-1  

RETURN 

END 


SUBROUTINE  PICKC  { NP T, Y,  PMAX ,PM I N, L PE AK , MTROG) 

TO  PICK  UP  PEAKS  AND  TROUGHS  FROM  A PROCESS 
Y — INPUT  PROCESS  FROM  WHICH  PEAKS  AND  TROUGHS  ARE  TO'  BE 
FOUND  __  _ . 

NPT  - NO  OF  INPUT  POINTS  OF  Y 

PMAX  — ARRAY  USE  TO  HOLD  SUCESSIVE  PEAKS  VALUE 
PMIN  — ARRAY  USE  TO  HOLD  SUCESSIVE  TROUGH  VALUES 
L PEAK  — NUMBER  OF  PEAKS  FOUND 
MTROG  — NUMBER  CF  TROUGHS  FCUND 

DIME  NS  ION  Y ( 1 f,  PMAX  I li,  PM  IN  I 1 > 

1=1 
L =1 
M = 1 

100  IF ( ( I+l ) .GE.NPT ) GO  TO  30 
25  I F ( Y ( I ) .GT. Y ( I+l  ) ) GO  TO  13 
20  I F C Y ( I J.GT.YC 1*1) ) GO  TO  15 
35  I F { I.GE.NPTI  GO  TO  100 

14  1=1+1 
GO  TO  20 

15  PMAX (L ) =Y( I) 

L=L  + 1 
GO  TO  100 

13  I FI  Yf  II.LE.YU  + 1 ) ) GO  TO  16 
33  IF( I.GE.NPTI  GO  TO  100 
18  1=1+1 
GO  TO  13 

16  PMIN(M)=Y( I ) 

M=M  + 1 
GO  TO  100 

30  LPEAK=L-1 
MTR0G=M-1 
RETURN 

end 
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SUBROUTINE  AREA  ( S, NPT, FU,F L J 
DIMENSION  SCI) 

COMMON  ANORML 

ANORML  — AREA  OF  THE  SPECTRUM 

FU  — UPPER  CUTOFF  FREQUENCY  IN  HERTZ 
FL  — LOWER  » * . « «• 

ANORML  = l .0 
PI2=2.0*3. 141593 
WUP=P I 2*FU 
WLR=PI2*FL 
DW=WUP/F LOAT ( NPT) 

NUMP=INT (WLR/DW) 

IF(NUMP.EQ.O)  NUMP= 1 
DO  45  1=1, NPT 
45  S(I)=0. 

DO  48  I=NUMP,NPT 
WI=FLOAT(I)*DW 
CALL  FUNSPT  (SP,WI) 

48  SII)=SP 

USE  SIMP SN  RULE  TO  COMPUTE  THE  AREA  UNDER  THE  CURVE 
NPP=NPT 

IF(NPP/2*2.EQ.NPP)  NPP=NPP-1 
CALL  SIMPSN  (S,NPP,DW»SAREA) 

ANORML=S  AREA 

WRITE (6, 60)  SAREA,FU,FL,WUP,WLR,(Sl 11,1=1,1281 
60  F0RMAT(10X,'AREA=',E14.T,5X,4C4X,F6.4)/8C  2X.F8.4J ) 
RETURN 
END 


SUBROUTINE  SIMPSN  ( AF , NPT, DSTEP, SI MP J 
DIMENSION  AFC  1) __ 

NPT  — NO,  CF  POINTS  FOR  INTEGRATION  (MUST  BE  ODD  NT, ) 
SIMP  — RE SULTI NG  AREA  OF  I N IE GRAT1  ON 

NP=(NPT-l)/2 

...  MN=NP-l  - 

0DD=0. 

EVEN=0, 

EMD=AF(1  J+AFJNPT  L - 

DO  10  1=1, NP 
10  ODD=ODD*‘AF  ( 2*  1 1 

00  20  I = 1,MN  - 

20  EVEN=EVEN+AF<2*I+1> 

SIMP=(4,0*ODD+2.0*EVEN+END)*DSTEP/3.0 

RETURN  _ 

END 
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SUBROUTINE  CROSS  ( ST OEV , PMA X ,PM I N, ILL , MMM , V ZERO , R AT  ID  I 

DIMENSION  COUNT  ( 51 ) , PMAX ( 1 ) » PM  IN ( 1 ) ...  

INTEGER  COUNT, SIGN 

TO  FIND  CROSSINGS  OF  A PROCESS  (INCLUDES  +VE  AND  -VE  SLOPS) 
WE  START  CROSSINGS  FROM  THE  FIRST  TROUGH 

STDEV  — STANDARD  DEVIATION  OF  THE  PROCESS  

PMAX  --  PEAK  VALUE  ARRAY  OF  THE  PROCESS 
PMIN  — TROUGH  VALUE  ARRAY  OF  THE  PROCESS 

LLL  — NO  OF  . PEAKS  ..  

MMM  — NO  OF  TROUGHS 
VZERO  — ZERO  CROSSINGS 

RATIO  — SET  TO  1 IF  CROSSINGS  FOR  THE  WHOLE  PERIOD  OF  THE 

PROCESS  IS  TO  BE  FOUND 

SET  TO  3600/PERI0D  IF  CROSSINGS/HOUR  IS  TO  BE 

FOUND....  . 

KL=1 

STDEV5  = STDEV/5.0 

DO  4 1*1,51  . _ 

4 COUNT ( I ) =0 
100  XXX=PMIN(KL) /ST0EV5 
YYY=PMAX (KL ) /STDEV5 
1 1 =X XX 
11=11+26 

IF(PMIN(KL).LE.O. ) GO  TO  11 
15  11=11+1 

11  KK=YYY 

KK=KK  + 26  ...  

IF{PMAX(KL J.GT.O  ) GO  TO  14 

12  KK=KK-1 

14  IFIKK.LE.5i  ) GO  TO  28  

33  KK=5 1 

28  IF ( II)  31,31,22 
31  11=1 

22  IF ( I l.GT.51 ) GO  TO  38 

20  COUNT! 1 1 ) = COUNT ( II ) + l 
38  11=11+1 

I F ( I I.LE.KK)  GO  TO  20 
25  XXX=PMIN(KL+1 )/STDEV5 

YYY=PMAX (KL ) /STDEV5  

IJ=XXX 
I J= I U+2  6 

IF(PM IN(KL+1 I.LE.0.0  ) GO  TO  21 
45  IJ=IJ+1 

21  K J=YYY 

K J=K  J + 26 

IF(PMAX(KL ).GT.0.0)  GO  TO  75 

74  K J=K J— 1 

75  IF (K  J- 51 ) 65,65,70 
70  K J=51 

65  I F ( I J ) 60,60,90 
60  IJ=1 

90  I F ( I J.GT.51 ) GO  TO  8 2 
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_80  COUNT! I J }=C  OUNT  ( I J)  + 1 
82  IJ=IJ+1 

IFUJ.LE.KJ)  GO  TO  80 
72  KL=KL+l . 

IF ( KL— MMM)  100,400,400 
400  WRITE! 6,500 ) 

500  FORMAT ( 1 HI ,5X NUMBER  OF  TOTAL  CROSSINGS  AT  DIFFERENT  LEVEL  OF* 
1,*  SIGMA,*  /6X, • INCLUDES  CROSSINGS  WITH  BOTH  *VE  AND  -VE  SLOP*// 
2 6X, 'SIMULATED  RESULT*//) 

177  Zl=-5.0  __ 

DO  66  1=1,51 

COUNT ( I J = I FI X( FLOAT (COUNT(I) )*RAT 10 ) 

WRITE! 6, 600)  Z 1 , COUNT !I) 

Zl=Zl+.2 
66  CONTINUE 
VZER0=C0UNT126)  _ 

600  FORMAT!  5X,  F4.1, ’ SIGMA* ' , 16, 2X, ‘CROSSINGS  • ) 

RETURN 
END  


C SUBROUTINE  SIMULT 

C PURPOSE 

C PERFORMS  SIMULATION  OF  RANDOM  SAMPLE  FUNCTION  BY  THE  USE 

C __  OF  FAST  FOURIER  TRANSFORM  METHOD  

C 

C USAGE 

C CALL  SIMULT  ! A »N  »NPT , NPT2, NPT4, WU, WL , I X, Ml , NSIGN , S , INV) 

C 

C DESCRIPTION  OF  PARAMETERS 

C A — A IS  THE  OUTPUT  SAMPLE  FUNCTION  IN  TIME  DOMAIN _ 

c IT  IS  A ONE-DIMENSION  ARRAY 

C N — N IS  THE  NUMBER  OF  SUBDIVISIONS  OF  A GIVEN  SPECTRUM 

C __ SS I W) » AND  N MUST  BE  A NUMBER  OF  POWER  OF  2 

C THAT  IS  N=2**M  ,M  IS  A POSITIVE  INTEGER 

C NPT  — NPT  IS  THE  NUMBER  OF  POINTS  OF  THE  OUTPUT  SAMPLE 

C __  FUNCTION  SPECIFIED  AS  INPUT  NPT=2**M1  WITH  

C 2*N  EQUAL  TO  OR  LESS  THAN  NPT 

C NPT2  — NPT2  IS  AN  INTEGER  INPUT  EQUAL  TO  2*NPT 

C NPT4  — NPT4  IS  EQUAL  TO  ONE-FOURTH  OF  NPT  ANp  MUST  ...... 

C BE  SPECIFIED  AS  INPUT 

C FU  — FU  IS  THE  UPPER  BOUND  FREQUENCY  IN  HERTZ  OF  THE 

C GIVEN  SPECTRUM  S$(W)  ,MUST  BE  SPECIFIED 

C FL  — FL  IS  THE  LOWER  BOUND  FREQUENCY  IN  HERTZ  OF  THE 

C GIVEN  SPECTRUM  SS!W)  ,MUST  BE  SPECIFIED 

C IX  —AN  ODD  INTEGER  INPUT  LESS  THAN  9 DIGITS, BY  CHANGING 

C IX, WE  WILL  GET  A DIFFERENT  SAMPLE  FUNCTION  FOR  A 

C GIVEN  N 

C Ml  — AN  POSITIVE  INTEGER  INPUT  SUCH  THAT  NPT=2**M1 
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C NSIGN  — AN  OPTION  PARAMETER 

C 0 USE  WHITE  NOISE  SPECTRUM  PROVIDE  BY  THE 

C SUBROUTINE  SIMULT  WITH  FREQUENCY  DOMAIN  WU  TO 

r WL  TO  SIMULATE  THE  SAMPLE  FUNCTION, AND  THE 

C SAMPLE  FUNCTION  IS  NORMALIZED  WITH  ME AN=0  AND 

C STANDARD  DEVIATION=l 

C 1 USER  WILL  HAVE  TO  PROVIDE  SPECTRUM  SS< W>  FOR 

C SUBROUTINE  FUNSPT  AND  THE  SAMPLE  FUNCTION 

C SIMULATED  IS  NOT  NORMALIZED  BUT  THE  MEAN  IS 

C STILL  ZERO 

C 

r REMARKS  

C THIS  SUBROUTINE  WILL  GIVE  A UP  TO  2**M1  POINTS 

C WITH  Ml  MUST  NOT  BE  LESS  THAN  3 OR  GREATER  THAN  20 

C 

C SUBROUTINES  AND  FUNCTIONS  REQUIRED 

C HARM  RANDU, FUNSPT 

c 

C METHOD 

C N 

C F ( T ) = SUM  SQRT (2.0*SS(W(K) ) *EXP ( I *PH II K )) *EXP ( I*W ( K ) *T ) 

C K=1 

C 

r 1=  IMAGINARY  NUMBER  T = T IME  PHI  =RANDOM  PHASE  ANGLE 

C AND 


C A(T)=SQRT(JWU-WU/N>*REAL(F(J)  ) 

C 

C FOR  N SUBDIVISIONS  OF  SPECTRUM,  FFT  ONLY  GIVE  N POINTS  OF 

C TIME  HISTORY,  THAT  IS  CNE  COMPLETE  PER IOD, THEREFORE  THE 

C RESOLUTION  OF  THE  PROCESS  IS  BAD.  TO  IMPROVE  THE  RESOLUTION 

C.  WE  HAVE  TO  ADD  SS(W(K)l=0  K=N+2 ,NPT  AFTER  SSIWUI 

C THE  TIME  INCREMENT  DTIME  1$ 

C DT  IME  = 1 • 0/W ( NPT ) W=FREQUENCY 

C NOTE:  SSCW)  HAS  PRESCRIBED  VALUE  ONLY  UP  TO  N POINTS 

C I F PHI  ARE  FIXED  UP  TO  K=N , THE  SAMPLE  FUNCTION 

C WILL  BE  THE  SAME,  THE  ONLY  THING  CHANGES  IS  DTIME 

C THEREFORE  WHEN  WE  INCREASE  NPT  WE  DECREASE  DTIME. 

C ..  A,  S, AND  INV  HAVE  TO  BE  DIMENSIONED  IN  TME  MAIN 

C PROGRAM  ACCORDINGLY. 

C 

C *****************************  *******************  *************_****** 


C 

SUBROUTINE  S IMUL T< A , N ,NPT ,FU ,F L , I X , Ml , NS IGN , I NV, S ) 
DIMENSION  AI1),M<3),  INI  V ( 1 » , S ( 1 > 

PI=3. 141593 
NL=256 
C 

C CONVERT  HERTZS  INTO  RADIANS 
C 

WU=FU*2.0*PI 
WL=FL*2.0*PI 
DW=WU/FLOAT ( N J 
NUMP  = INT  C WL  ADW) 

WL=DW*NUMP 


-46- 


n o o | r»  jr»  r»  | r>  o ;r>  j 1 I I j i r>n|o 


IF(NUMP.EQ.O)  NUMP=1 

SQR00T=SQRT(2.0*WUZ(  ( WU-WL1  *FLOAT  (N  ))  I 

Ml  1) -Ml 
M<2)=0 

M ( 3 ) =0  _ _ . _ : 

if$et=i 

00  5 1=1 »NUMP 

A(  1*2-1 1=0.0 

A l 1*2) =0*0 
5 CONTINUE 

NN-N+2  _ 

00  10  I-NNtNPT 
A l 2* 1-1 )=0.0 

_ A(  2*1 > = 0.  0 _ 

10  CONTINUE 

IF(NSIGN.EQ.l)  GO  TO  40 

TO  GENERATE  RANDOM  PHASE  ANGLES  PHI 

DO  30  J=NUMP1N  

CALL  RANDU  (IX, IY,YFL) 

IX=IY 

P H 1=2  «_0*  P I * YF  L _ __ 

A( 2* I+1)=C0S ( PH  I ) 

A( 2* 1+2 ) = S INI  PHI  ) 

30  CONT  INUE  _ _ _ 

GO  TO  60  

_40  00_50_J=NUMPjiR __7. ' 

W1=FL0AT(J1*DW 
CALL  FUNSPT  (SP,W1> 

_ CALL  RANDU  K XX*  IY.YELI  _ _ 

1 X=I  Y 

PHI=2.0*PI*YFL 

_ A(2*J±lL=SQR_I(2iO*SP}*COSJPHIl __  

A<  2*J*2l  = SQRTI2.0*SP)*SINlPHn 
50  CONTINUE 

START  S I MU LAtT N G SAM P L E FUNCTION  FOR  A GI VEN  SET  OF  PHI 

60  CALL  HARM  (A»M»  INV, S,  I FSET«  IFERR  1 

IFlNSIGN.EQ.il  GO  TO  80 

SAMPLE  FUNCTION  FROM  SSCW)  ECUAL  TO  WHITE  NOISE 

DO  70  1=1, NPT 

_AI_n=A(2*I^l_l*SQR00T _ _ 

70  CONTINUE 
GO  TO  100 


SAMPLE  FUNCTION  FROM  A SPECIFIED  SSIW) 

80  DO  90  1*1, NPT 

A(  I l=A  <2*1-1 l*SQRT< DW) 

90  CONTINUE 

100  RETURN  _ 

END 
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SUBROUTINE  fUNSPT  (SP,W1) 

C . . . ..  ..  . ........  . . .. 

C SS(W)=...  IS  A FUNCTION  STATEMENT  NONEXECUTABLE  AND  MUST  BE 

C SUPPLIED  BY  THE  USER.  SS(W)  IS  THE  SPECTRAL  DENSITY  FUNCTION 

C 

COMMON  A NORM L 

SS(W)=P12*PHIO/(1.0+(PI2*W/F8)**A)  /ANORML 
P 1 2=1 •/ ( 2. *3. 1A1 593 ) 

PH!0=30 
FB=0 .031 
A = 2.  6 
SP=SS(W1 ) 

RETURN 

END 


C 

C 

C 

C 

C 

c 

c 

c 

c 


SUBROUTINE  PLOT  (MPT  »H *FUN, F ACT, NSTEP ) 

NPT NO.  OF  POINTS  TO  BE  PLOTTED 

FUN FUNCTION  TO  BE  PLOTTED 

H TIME  INTERVAL  (TO  2 DECIMAL  PLACES) 

REMARK:  TITLE  AND  DIMENSION  OF  FUNCTION  HAVE  TO  BE  PRINTED  OUT  BEFORE 
CALLING  PLOT 

REMARK:  IF ' JJCT-O..THEN  SCALE  FACTOR  IS  COMPUTED  BASED  ON  THE  GIVEN 
FUNCTION, OTHERWISE  SCALE  FACT QR=  FACT 


01  MENS  I ON  ALINE  ( 120)  , FUN  (NPT  ) , DA  S ( 1 20  ),,  SCAL  E ( 1 1 ) 
DATA  DASH/,-,/,STAR/,*,/, BAR /' I'/, BLANK/*  •/ 

C 

C ..  DETERMINE THE  SCALE  FACTOR 

IFIFACT.EQ.O. ) GO  TO  30 
FACTOR=FACT 
GO  TO  80 
30  AMAX=1. 0E-70 

DO  M 1 = 1, NPT, NSTEP 

l F ( ABS( FUN( I ) ).GT.AMAX ) AMA X=ABS ( FUN  ( I) ) 

AO  CONTINUE 

FAC  T OR  = A M A X / 5 0, 0 


C START  TO  PLOT 

C 

80  DO  95  1=1,11 

95  SCALE!  I )=-60.+l  0.*F1.0AT(  I ) 

DO  100  1=1,120 
100  OAS(  I)  = DASH 

DO  105  1=15,115,10 
. 105  DAS(  I ) = BAR 

DO  110  1=10,120 
110  AL  IN E ( I ) = BL ANK 
AL  INE( 65 ) =BAP 
WRITE (6, 200 ) FACTOR 
WR ITE  ( 6, 302  ) SCALE 
WR I TE ( 6, 303 ) (DAS( I ), 1 = 10,120) 


-48- 


o o io  r>  o r>  o r»  io 


DO  120  1=1 ,NPT,NSTEP 
T IME=H*FL0AT*I-1 ) 

J*P=FUN*  I )/FACT0R±65.  

ALINE*NP)=STAR 

WRITE (6, 304)  T IME , ( AL INE I J) ,J=10,120) 

ALINE*NP)=BLANK  __  _ 

IFINP.EQ.65)  AL  INE*  65 ) =BAR 
120  CONTINUE 

200  FORMAT*///,  50X,  ‘SCALE  FACTOR**,  _EH.3) 

302  FORMAT  < 1 HO , 2X , * T IME • ,5 X ,5  * F4.0 , 6X ) , IX , 6 < F3. 0 , 7X ) ) 

303  FORMAT  ( 1H  ,8X,114A1) 

304  FORMAT  *JH  ,F8.4,111A1)  _ 

RETURN 

END 


SUBROUTINE  P ICKMX* A , NP ,PE AK , MPK,TROUG , NPK ) 

PICK  UP  THE  DIFFERENCES  BETWEEN  SUCESSIVE  MAXIMUM  RISE  AND 
MINIMUM  FALL  ( 2ND  METHOD  OF  PICKING) 

A — INPUT  ARRAY  FROM  WHICH  RISE  AND  FALL  ARE  TO  BE  P ICKED 

NP  — NUMBER  OF  POINTS  OF  INPUT  ARRAY  A 

PEAK  — _ ARRAY  USE  TO  HOLD  THE  RISE  AND  FALL  DIFFERENCES  . .. 

MPK  — NUMBER  OF  RISE  AND  FALL  DIFFERNENCE  FOUND 

_ DIMENSION  All)  ,PEAK(1)  , TROUGH) 

AMAX=C. 

AMIN=0. 

IFLAG-1  _ 

IF(A(l).LT.O.)  I FLAG=0 
I PK=1 

ITK=1  _ 

NST*  1 

IF* IFLAG.EQ.O)  GO  TO  20 

DO  10  1 = 1 ,NP  

IF C A*  D.LT.O.O)  GO  TO  15 
10  CONTINUE 

_15_NST=I  

iflag=o 

20  DO  100  I =NST ,NP 
__  IF*  IFLAG.EQ.O)  GO  TO  40  __ 

IF* A*  D.GT.AMAX)  AMAX=A*  I) 

IF* A 5 IH  ).LT.O.  ) GO  TO  50 
GO  TO  100 
50  I FLAG* 0 

PEAK* IPK)=AMAX 

IPK= IPK+1 

AMIN=1.0 
GO  TO  100 
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non 


40  IF(A(  D.LT.AMIN)  AM  I N=  A(  I) 
11=1  A I mj.GT.O.  ) GO  TO  70 
GO  TO  100 
70  I FLAG=1 

TR0UG< ITKI=AMIN 
I TK= ITK  + 1 
AMAX=-1,0 
100  CONTINUE 
MPK= IPK-1 
NPK= ITK-1 
RETURN 
END 


SUBROUTINE  STATIS  ( A , NP , STD , RMS ) 

DIMENSION  A ( 1 » 

TO  FIND  THE  STANDARD  DEVIATION  AND  RMS  VALUES  OF  A PROCESS 


AS=0. 

AS2=0. 

DO  10  1=1,  NP 
AS=AS+A(  I ) 
AS2=AS2+A» I )*A{  I ) 

10  CONTINUE 

A S= AS/FLOAT (NP) 

AS 2= AS 2/ FLOAT (NP) 
STD=SQRT(AS2-AS*AS) 
RMS=  SQRT ( AS  2 ) 

RETURN 

END 


SUBROUTINE RANDU (IX, IY, YFL I 

IY=IX*65539 
IF(IY)  5,6,6 

5 . IY-IY+2147433647+1 

6 YFL=IY 

YFL=YFL*0. 4656613 0-9 

RETURN - - 

END 
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SUBROUTINE  H AR  M ( A ,F  , IN  V , St  IFS ET  , I F ER R ,NPT 2 ,NP T 4 ) 

DIMENSION  A (NPT 2 ) , I NV ( N PT4 > , S (NPT4  J , N ( 3 ) , H 1 2 ) t NP ( 3 ) ,W { 2 ) , W2 ( 2 ) , 
IW3 ( 2 ) 

EQUIVALENCE  ( M , N 11 ) ) , CN2, N { 2 ) ) , ( N 2, N * 2 ) } 

10  I F ( IABS(IFSET)  - 1)  S0C,900,12 

12  MTT=MAX0(M(1) * M (2  ) * M(  3 ) ) -2 
RP0T2  = SQRT( 2 . ) 

IF  (NTT-MT  ) 14,14,13 

13  I FERR= 1 
RETURN 

14  IFERR=0 
M1=M  1 1 ) 

M 2=  M ( 2 ) 

M 3=  V ( 3 ) 

N 1=  2**M 1 
N2=2**M2 
N3=2**M3 

16  IFUFSET)  18,18,20 

18  NX=  N1*N2*N3 
FN  = NX 

DO  15  I = 1,NX 
A ( 2*  I—  1 ) = A « 2 I — 1 ) / F N 

19  A (2*1  ) = -A  (2*  I )/FN 

20  NP ( 1 )=N 1* 2 

NP  ( 2 ) = NP ( 1 )*N2 
NP(3)=NP(2)*N3 
DO  250  I D = 1 ,3 
IL  = NP  C 3 J— NP (ID) 

I LI  = IL  + 1 

mi  = y ( i d ) 

IF  (y I >250,250,30 
30  I DI F = NP ( I C ) 

K B I T = N P ( IC) 

MEV  = 2*(yi/2) 

IF  (M  - PEV  >60,60,40 
40  KBIT=KBlT/2 
KL=KB  IT-2 

DC  50  1=1, 1 LI, I OIF 
KLAST  =KL+  I 
DO  50  K = I , KLAST, 2 
KD=K  + KB  IT 
T=A (KC) 

A(KD>=A(K  )-T 
A(K)=A(K)4T 
T=A (KC  + 1 ) 

A(KC+1)=A(K+1)-1 
50  A(K+1)=A( K+l J4T 

IF  (MI  - 1)  250,250,  52 
52  LFIRST  =3 
JLAS7=1 
GO  TO  TO 
60  LFIRST  = 2 
JLAST=0 
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70  DO  240  L=LFIRST ,R I ,2 
J JD I F = K E I T 
KBIT=KBIT/4 
KL=KB IT-2 

DO  80  1=  1 , IL1 , IC  IF 

KLAST=I +KL 

DC  80  K=I,KLAST,2 

K 1=K+  KB  IT 

K2  = K1+KE  IT 

K3=  K2  +K E I T 

T =A  ( K 2 } 

A (K  2 ) =A  { K )-T 
A ( K ) = A ( K ) 4T 
T =A ( K 2+ 1 ) 

A(K  2+  1 ) = A (K  + l )-T 
A (K+II=A (K+l ) +T 
T = A (K3 ) 

A(K3)=A(K1)-T 
A ( K 1 ) = A ( K 1 ) +T 
T=  A (K3  + 1 ) 

A(K3+1  )=A(K1  + 1 )-T 
A(K1+1)=A(KH-I)4T 
T=A ( K 1 ) 

A ( K 1 ) = A ( K )-T 
A ( K ) = A ( K ) +T 
T = A ( K 1 + 1 ) 

A ( K 1 * 1 J = A (K  + l )-T 
A(K+1)=A(K+I)+T 
R=- A ( K 3+  1 ) 

T = A (K3  ) 

A ( K 3 ) =A  ( K 2 )— R 
A(K2)  = A(K2)+R 
A ( K 3+ 1 ) = A {K2+1 J-T 
80  A(K2+1)-A(K2+1)+T 

IF  (JLAST)  235,235,82 
82  J J= J JDI F +1 
I L A ST=  IL  + J J 
DO  £5  I = JJ, I LAST, IDI F 
KLAST  = KL4I 
DO  85  K= I , KLAST  ,2 
K 1 = K4KBIT 
K 2 = Kl  + KBIT 
K 3 = K24KPIT 
R =-A ( K 2 + 1 ) 

T = A ( K2  ) 

A ( K2 ) = ACKI-R 
A ( K ) = A(K)4R 

A ( K 2+ 1 ) = A (K4IJ-T 
A ( K 4 1 ) = A { K + I ) 4T 
AtoR=AlKl)-A(Kl+l) 

AWI  = A ( K 141 J 4 A ( K 1 ) 

R=-A I K3 )-A ( K34  1 } 

T A ( K 3 ) - A (K341 1 
A(K3)=(AWP-R)/RCCT2 
A (K  34  1 ) = ( AWI-T  I/RC0T2 
A(K1)  = ( AV»P4R)/RC0T2 
A (K141 ) = ( AWI4T  J/RC0T2 
T=  A ( K 1 ) 


A (K 1 )=£ (K )— T 
A (K ) =A { K ) +T 
T=A ( K 1 + 1 ) 
A(K1+1)=A(K+1)-T 
A (K+ 1 ) = M K + l ) +T 
R=—  ACK3+1) 

T=A (K  31 
A(K3)=« (K2)-R 
A ( K2 ) =A { K2 ) +R 
A(K3+1) = * (K2+1J-T 
85  A(K2+1)=A (K2+1J+T 

I F { JL £ ST~ 1 ) 235,235,90 
90  JJ=  JJ  ♦ JJOIF 
DO  230  J = 2 » JLAST 
96  I=INV(J41) 

98  I C=  NT—  I 

wm=suc) 
w ( 2 )=sm 
12=2*1 
I 2C=NT-  12 

IFU2C)  120,110 ,100 
100  W2C 1 1 = SC  I2C  ) 

H2 ( 2 ) =S ( 1 2 ) 

GO  TO  130 
110  W2 ( 1 ) = 0 • 

W2(  21=1. 

GO  TC  130 
120  I 2CC  = I2C+NT 
1 2C.  =—  1 2 C 
W 2 ( 1 )=-S ( I2C) 

W2  ( 2 1 =S  ( 1 2CC  ) 

130  13=1  + 12 
I 3C=NT-I3 

IFU3C)  160,150,140 
140  W3( 1 ) = S ( I3C ) 

W3( 2 ) =S(  13) 

GO  TO  200 
1 50  W3(  11  = 0. 

W3(2)=l. 

GO  TO  200 
160  I 3CC=  13C+NT 

I F { 13  CC J 190,180 ,170 
170  I 3C=-  I3C 

W3  ( 1 ) =-S  ( 130 
W3( 2 ) =S ( 13CC ) 

GO  TC  200 
ieo  W3C l ) =- i . 

W3C  21  = 0. 

GO  TC  20C 
190  I 3CCC=N  T+  I3CC 
I 3CC  = -I3CC 
W31  1 )=-S  C I3CCC ) 

W3 ( 2 ) = -S ( I3CC ) 

200  I LA  ST=I L+ J J 


DO  220  I=JJ,I LAST, IDIF 

KL * 57  =KL+  I 

DO  220  K=I,KLAST,2 

K 1=K  + KB IT 

K 2 = K 1 + K B I T 

K3=K2*KB IT 

R=A (K2)*k2( 1)-A (K  2+ 1 } *W2 ( 2 ) 
T = A(K2)*W2(2)-*A(K2  + l)*lft2(l) 

A (K  2 ) =A  ( K )-R 

A(K)=A(K)+R 

A {K2+ 1 ) =A (K+l )-T 

A(K+1)=A(K+I)+T 

R=A(K3)*l»3(l)-A(K3+l)*h3(2l 

T=A(K2)*fc3(2MA(K34l)n2(l  ) 

AWR=AlKl)*Wll)-A(Kl-tl.)*V\(2) 

AWI=A<K1)*W(2)+A{K1  + 1)*M1) 

A ( K3 ) = A hR-R 

A { K3+ 1 ) =A WI -T 

A (K1)=AV»R*R 

A IK  1+  1)  =AViI  +T 

T=A(K1) 

ACK1 )=A(K  )-T 
A(K)  = A(KMT 
T=A(K1  + 1 ) 

A(K141)=A(K+1 )— T 
A (K+ 1 ) -A ( K+l) +T 
R=-A  ( K3+  1 ) 

T =A  C K 3 ) 

A (K3  ) = A(K2)-R 
A (K 2 ) =A { K 2 ) +R 
A (K  3 + 1 ) = A (K2+1J-T 
220  A(K2+ 1 )=A (K2+1 )+T 
230  JJ=  J JC  I F 4JJ 
235  JLAST=44JLAST+3 
240  CONTINUE 
250  CONTINUE 
NTSC=NT*NT 
M3MT=F3-y T 

350  I F ( M 3RT  ) 270,360,360 
360  I G03= 1 

N3VNT=N  2/M 
MINN3=NT 
GO  TO  380 
370  I G03=2 
N3VNT= 1 
NTVN3=NT/N3 
MINN2-N3 

380  JJD3  = NTSQ/N3 
M2MT=M2-MT 

450  IF  (M2MT )470,4fcC,460 
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460  IG02= 1 

N2VM=N2/M 
MI NN2=NT 
GO  TC  480 
470  I GO  2 = 2 
N2VNT  = I 
NTVN2=N7/N2 
M INN2=N2 

480  J JD2=NT  SC/N2 
M 1MT=M 1-MT 

550  I F ( M1MT ) 570 ,56C .560 
560  I GO  1= 1 

N1VNT=N1/NT 
MINN1=NT 
GO  TC  58C 
570  I GO  1=  2 
N1VNT=1 
NTVN 1 =NT / M 
M INN1 =N 1 

580  J JD 1=NTSG/N1 
600  J J3= 1 
J=1 

DC  880  JPF3=1 ,N3VNT 
I PP3= INV ( J J3 ) 

DO  870  JP  3= I » M t NN3 
GO  TC  (610,620) ,IG03 
610  I P3= I NV (JP3)*N3VMT 
GO  TC  630 

620  IP3=INV(JP3)/NTVN3 
630  13= ( 1FP3+IP3) *N2 
700  J J2= 1 

DO  870  JPP2=1 , N 2V  NT 
IPP2=INV( JJ2J+I3 
DC  860  JP2=1,MINN2 
GO  TC  ( 710,720) ,IGC2 
710  IP2=INV(JP2)*N2VNT 
GO  TC  730 

720  IP2=INV( JF2 1/NTVN2 
730  I 2= ( I PP 2+ I P2) *N1 
800  JJ1=1 

DO  860  JPP1=1 ,N1VNT 
IPP1=  INV  C JJD  + I2 
DO  85C  JP  1= 1 ,M INN  1 
GO  TO  (£10,820), I GO l 
810  IP1=INV( JP1)*MVNT 
GO  TO  830 

820  IP1=INV(JF1)/NTVN1 
830  I = 2 * ( IPP1  + IPD  + 1 

IF  (J-I)  £40,  £50,  £50 
840  T=A ( I ) 

A ( I ) = A ( J ) 

A ( J ) =T 
T=A (1*1) 

A ( I + 1 1 = A ( J+ 1 ) 

A ( J *1  ) = T 
850  J = J+  2 
860  JJ1=JJ1+JJD1 
870  J J2=J J2+J JD2 
880  JJ3  = JJ3+JJD3 
850  IF(IF£ET)891,895,£95 


a <31  00  892  I = 1 f NX 
8 9 2 A ( 2 * I ) = -A(2*I ) 

895  RETURN 

90C  M T=  MA  X0  ( M ( 1 ) , M ( 2 ) , M ( 3 ) ) -2 
MT  = MAX0(2»MT) 

904  IF  (MT-18)  906,906,13 
906  I FERR=0 
NT=2**MT 
NTV2=NT/2 

910  THETA=. 7853981624 
J ST  E P=N  T 
J 01  F=NTV2 

S ( J D I F ) = S INITFETA  ) 

DO  950  1 = 2 , MT 
THE  T A = THE  T A/2 • 

JSTEF2=JSTEP 

JSTEP=JDIF 

JDIF=JSTEP/2 

S(JDIF)=SIN(THETA) 

JC1=NT-JDIF 
S( JC1)=C0S(THETA) 
JLAST=NT-JSTEP2 
I F ( Jl A ST  - JSTEP)  950,920,920 
920  DO  940  J = JSTEP  ,JLAST, JSTEP 
JC=NT-J 
JD= J+JDIF 

940  S( JD)=S( J)*S( JC1)+S (JDIF )*S (JC) 
950  CGNTINUE 
960  MTL  EX  P=N  T V2 
LM1EXP=1 
I NV ( 1 ) =0 
DO  980  L = 1 » MT 
I NV (LN1EXF  + 1)  = MTL  EXP 
DO  970  J=2,LM1EXP 
JJ=J+LM1EXP 

970  INVl J J)=INV(J)+MTLEXP 
MTLEXP=MTLEXP/2 
980  LM1EXP=LM1EXP*2 
982  IF  C IFSET  ) 12  ,895 ,12 
END 
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APPENDIX  B 


//  E<EC  FQRTGCG 

//F5RT.SYSIN  DD  * , DC  B=R  LKSI  7.  E=3  V+  0 
C ARFL1 

C SENSITIVITY  STJDY  15  AVAILABLE 

X-Q  =— XA  -K  .X4  RE-SL1LD  X-S-  AVAIL  ABLE. 

C THREF  LOADS  ( COMPOS  1 1 OR  SINGLE  GAUSSIAN) 
DIMENSION  A A (10)  , AAAI  10) , BAB  (10  ) 
DIMENSION  A ( 1501)  ,R  (15  01  ) ,H(  15  01  ) , 


ARE  AVAILA3LE 

T H ( 1 5 0 1 1,(1(1501) 


-aiME-NS-inN--P.P4a| ,SGC(3  U-RR(  3), &Mf-TA(.3 1 

DI MENS  ION  DUMY(1 5 01 ) , TP ( 5 0 ) ,P ( 50 ) , PS ( 50 ) , RK ( 50 ) , Cl J ( 50) , 31 J( 50  ) 
DIMENSION  ;MP(3),SGS(3I,VX(3),AN(3),UP;(3)».S8(3) 

DIMENSION  AAAA(IO) 

4L0ALT4-N-U&- 


£ *?  #'  5*: rir-'*: $. tfr  <r  ■& SI. ?£;  jfr £■  3*  $ *>  *%  :$t  < fr  4 i*  t # s? 

■:  *5.’  t*s  ,*•  ^ ^ * 1J5  *c  M F A ^ V A L J P S *"  # * * * < • * j: 

^*-*$*fcj'*  =RA;TURE  DAT  A **, 


Qf  S- 

-4=E 


f ■ if !-  $■  i:  t 15: 4;  S:  fc  f |c  |r 

(J:  ft  Js  fc  «jp; '/  >ifc  j>:  :f  sft  j$:  ift 

'St  J fc  'fc  t*  iV  rf.  *'<  ?*•'  $ !>'•  t 1*  f? 


=-2~5- 


BE=A.O 
DKTH  =1.5 
EFD= 10  .0 
C=3  . OE-7 


KG  = 60.  0 

£*•!*!*>;  !;  hi . js  (s  INITIAL  CRAC<  DATA  * * *r***«k*  * t('<8  fc  6 * £ * 

A 0=  0. 04 

ALPHA  =4.-0 

BET  A=30000.0 

EHkM1  (*!**$*  4ifc  (■:  * LOAD  DATA  (s*  * (•*  ******  ♦♦****♦**  ******* 
C CM 5 (1.0=C0MP3S IT  GAJSSUN,  0.0  = SINGLE  ONE) 

AMP  ( l ) =1  .0 - - - - ....  ' 

CMP ( 2)  =1.0 
CMP (31=1 .0 

C up;  (1.0=UPPER  CUT  SINGLF  GAUSSIAN,  3.0  = C3M>LETE  ONE) 

UP  C 114  =0  .0 __ 

UPC  (2)  =0.0 

up; ( 3) =0.0 

C E-  = UPPER  CUT  LEVEL  (EL  * X3  ) 

— FL=3-»0 

G= 10 . 0 
AA(4)=  115.0 
AAA  ( 4)  =115.  0 

A-.AAjU-U--BD.-0 

SAB (4) =3 .0 

C AN ( I ) NUMBER  OF  LOAD  CYCLES  PER  HOUR 
AN ( 1)  = 600 .0 

AN  (.2  1 = 600.0 

AN(  3)  = 60.0 
ANA=0. 5 
PP (1 ) =0. 4975 


c sg; 

PP(3)=0.5 
STANDARD  DEVIATION 

sg;(  n=o.o7*G 

RGC.  (?l=O.QBAr; 

OF  COMPOSITE  GAUSSIAN 

C SGS 

SGC  ( 3 ) =0.1*G 
STANDARD  DEVIATIONS 
SGS ( 1 ) =1 .5*0 
SGAL2J  =1  . 511G 

(SINGLE  GAUSSIAN) 

SGS  ( 3 ) =1  . 5 G 
X0=1 .O^G 
7R=1.5 

.......  . 
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r**ft  k**:**  M » MAT  E R I AL  DATA  *■  •'*■’•*  ■ * ■’ * t t' ^ N* **■>, >:  >?* 

V0= 0.056 
AMJ0  = 5.7*-G 
P 0 = 5 . 7N'G 

r **■*•■  ♦ •'*&■■! si  •£■  * FATI  SA  = = TATA  * it  • t i ft  < i frk  * * (■  *• t ,.  t i u.  y 

F5F=1.0 
ZFTA=0.43 
AS=  7 .0 

F ;%  * t:  & K'  >V  St  ft  ^ •'  !'■  ■’  it  rft  *■»:  K s’  F +■  »■  <c :'  *■*’(:'  t i ft’  e::  ■:  ’.'S  v ft  r\  3:;S-  t\*V  »v!*>  ft  { ft;  ft;  tc 

C TT  TOTaLFOGH T HOUR  : 

TT= 15000.0 

INSPECTION  DATA  tv**  **  fr**  * tj-.w**  * A * 

AI  -0.0?. _ _ ! 

A2=  0 .3 
RM=0.l 25 
U=l.  0 

■T.  NM  NDMBER  OF  INSPE C JJLflN S 

C TT  / ( MN 4- 1 ) MUST  8F  IMT  EGER 
C (TT/{NN*in/<TT/(NPT-in  MUST  BE  INTEGER 
NN  = 0 

fe  ;f  •&  i*r  t} : *•  i$  j%-  ■*!  j?1 it  ?V>?  * jf:  tS- 1-  $•  f.  + fr  £ 9 :&s»  ^ ^ i»-  i £ 

p i = 3 . 14159  * 

C NPT  BAS  IC.  0 ATA  POINTS 

C EACH  DIMENSION  SIZE  0 F ( A , R , H,  TH,  OUMY  , N )M  JS  T BE  G ,E  . NPT 

r.  ITT  /(  N 3 T- 1 I)  M u ST  RF  I NT  - PER 

NPT =15 01 

C NO  I DATA  POINTS  UNDER  NO  INSPECTION  CASE 

C EACH  DIMENSION  SIZE  DF(TP  P , P S » P K , C I J , 01  J > M JST  BE  G.E.  NDI  AND  (NNftll 

-C  (TT/NOn  MUST  BE  INTEGER 

C ( TT/NOII  /(  TT/(NPT-in  MUST  BF  INTEGER 
NO  I = 30 

C ; + #■«•  ft  M'fc  ft  * * * iftfe  ft  * FLEET  SIZE  t f *tft  *»*••<:*  < *’ 4- :*  *,  *;•  >f  <•** 


■ >'■  sj. ;»  j»  j' it  ift  !:  *■ 


£***  AtftjftlM  *******  ft*  A ft.lftft  ft  ft  * ft:>  ft  A k'ft  (t  V ■'  4 f Wt  * *,  ip!  FF?Ot  iFiETTf  ft  E i 
r *scT|t  *?  h a a **  **••?:*  A**:  *♦  ***«•«!  ift*  t * <t*4  * »'  t |c  ► j-.  k ft  k k It  -Mi  y it  kit 

C SUB  PROGRAM  FOR  SENSITIVITY  STJDY 


501  FORM  AT  (4  FID. 0) 

WR I TE ( 6,650)  XY , Y7 1 , YZ 2 , YZ 3 
6>  0 FORMAT (1H  ,///,’  X Y, YZ 1 , Y Z 2 , YZ 3 » , 4 E 1 5 . 7 ) 
IF  f XY.LE.O.O)  GO  TO  32  2 . 


IF C XY.EO  .1.01  FF  = YZ 1 
IF  ( X Y . EQ.2.0  ) DKT  H = YZ  1 
I = ( XY.EQ.3.  0)  Ec 0 =YZl 

IFIXY.FD .4.0 ) C= YZ  I 

If  (XY.E0.5.0I  AO  = Y Z1 
I F ( XY.EQ.6. 0)  AN0  = YZ l 
IF (XY. E0.7. 0 ) ANA=Y71 


1=  ( XY.  ED  • 9.  0)  V0  = YZ  1 
IF(XY. EQ.1D.0)  Z ET  A=  Y Z 1 
IF(XY.  EO.ll.O)  A S = YZ l 


IF(XY.NE.13 .0)  GO  TO  721 
PP(1)  = YZ  1 
PP(2)  = YZ2 

PP  ( 3 i =YZ  3 

721  CONTINUE 

IF(XY.NE.14.0>  GO  TO  722 
SGC.(1)  =YZ1 
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s GC  ( 2 ) =Y  Z2 
S3:  (3)  =YZ3 
7??  CDNIT  INUF 

IF(XY.NE.15 .0) 

...  G-YZ1 

SGC ( 1)=0.07*G 
S3C  ( 2 ) =0 . 1 8*  G 
SGC(  3>=0.l*G 

— 

R 0=5. 7* 3 
AM'J0  = 5 .7*G 
7?  3 CONTINUE 

I.F4.-XY-..NE^16..-Hi 

AMJO= Y Zl 
R0=YZ1 

7?  4 COMTIMUE 

XF4XY^FCUL7*0-i- 

IF( XY. EQ. 18. 0) 
IF(XY.E3.19.0) 
IF (XY. EQ.23 .0) 
IJrJJCV.  Efl*21*43J- 


GD  rn  723 


-Ga..-XQ-2212A- 


-A_-B-HA=-YLZ-i- 


8 E T A =Y  7 1 
F S F = YZ  1 
RE=YZ1 
■A A ( Al=Y_ZL 


4Z2S- 


IF(  XY.E3  .21  .OJ  A A A ( 4 ) = YZ  1 
IF ( XY. FQ.2? .0)  AA  A A ( 4 ) =Y  71 
T F ( XY.FQ.23. 0)  7R=YZ1 
IF4XY.NE.24-.  01  _G1_H_L25__ . 

sssm  =yzi 

SGS( 21=YZ2 
SGS  ( 3 ) =Y  Z3 
COHTiNllE 


IF(  XY .NE .25 .01  33  T3  7 26 
CMPm  =Y  7 1 
CMP<  2) =YZ2 

£Jl?_LALsYZA 

7?  6 CONTINUE 

I F ( XY.  EQ  . 26.  01  A1=YZ1 
IF  C XY.EQ.27. 0)  A2=YZ1 
IF4-Xy-.FQ.-23^01  .RM  = YZ4 


I F f XY.  EQ.29.  0) 


0) 

i 01 


72  7 


J1L_£J.=YL2I. 


r sj  3 

C***  r 


IF(XY. EQ.30 
IF(XY.NE.31 
U3C(  l)=YZl 
IJPC  ( 2 I = Y Z2 
up: < 3)  = YZ3 
C3MT  mjF 
-l  E4XY.ECL.  32 
PROGRAM  FOR 

4?  $ 4k  iiz  ok 

JKL  =0 

_X.BF=3E 

I P F = B E 

I F ( C M P { 1 ).F3  .1.0) 

I F ( CMP  a I.FQ.O.O) 
-1-4  : MP  ( 2. 1 ..Eil.1  .jQI 


J = YZ  1 
NM  = Y7.1 
GO  TO  727 


SENSI  T I VI  TY  STUDY 

&£  £ -■  h $ -|7  f ' $ iff  tr  £ S2’:Y';; 


*>■»•  END 

. s;  t <• 


: if. * * t £ Sr  S > i fc 


S3  ( 1>=AA<  IBEH3AB(  IBE)*SGC(  1 »*KBE 
SB(1)=AAU8E)  *SGSm*t<3E 

SB  ( 2 ) = A A A ( I BE  )4  B AB  < IBEHSGCI2  )-'-*KBE 


IF ( CMP (2  ) .Eg  .0  .0 ) S3 (2 | = AAA ( I BE) 
IF  { CMP  (3  ) .EQ  .1  .0)  S 9(3  ) = AA  AA  ( I RE)' 
IF(CMP{3).FQ.0.0) SR(3) =AAAA( IRE) 
ZB3  = (ZRiG)  : -;KRE 


» S3S(  2 )?-K3E 
BAB(IBEI*SGC(  3)**<3E 
* S G S ( 3 ) ! * K R F 


0 = AN(  l I*  PP(  1 )*  SBC  1 )+AN  12  )«PP(2  M SB  <2  M-AM(  3)*PP  (3)  *SB(3> 
1 f AM  A v7  BB 

AC  = 2 . / PI  * (K  C/R3  ) -’"2 
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DT  = TT/FLOAT(  NPT-l  ) 

C 

ANS  = PP  (1  KAN  ( 1 l+PP  ( 2 ) f AN  { 2 ) +P  P ( 3 ) "■  AN  ( 3 ) 

r ans=cy:lfs  pfr  hour 

C - S $ $ - THR  FA  HOI  n 4 0 1 RET  A At  I ) MD  I SAN!  _tl$ 

TF( EF. EQ.2.3  ) FF  = 2 . 3 I 
TF(  FFD.E'O.  2. 0)  FFD=2.0l 
C 

- A N T-H  .=  A L O G 1 3 (X)  *FF  tA  L 0 3- 1 0 ( DXXHi 

r 

CDST-ANTH-FFn?  AL3G10OKTH) 
r d s = 1 3 .n**r.n$T 

C : 

AH=  { OKI  H/  ( Q/ANS  )**{  1.0/3FI)  : ' 2'  <2. 0/P  II 

r 

T3  = 2.0/(  2.  O-FFO)  / CDS/ Or  UEFD/BF)  *ANS^fr  ( ( EFD-BE  )/BE  ) 

TR-TB/l?  1/2 .0)  : ( EFO/2.3) 

TR  = T B*  ( ATH%  is  ( (2.3  -F;FD  ) /? . 0 ) -A0«  * ( ( 2 .0- EFD  K 2 .3  ) ) 


C IN 


yn  43 


03  8 1=1  , N>T 

T-F.l.OA  T(  1-1  J-DT 

I F { TR.GT.O.  0)  GO 
CAFE  OF  ATH.LT.AO 
A ( I I =A0<!'  * ( ( 2 ,0-EF  ) /? 
_1U1=UU±U2 
1 » (PI/2.0) 


IF( A (I) 

A ( I )=A (I  »**{ 2 
GO  TO  43 


3 ) 

Q-XEl/Z,  0U:H3-*-..(E=/BE)  *ANSt  * ( (BE-EF)  / 3E)  » r 
t(  EF/2.0) 

LT. 0.00031 ) 33  TO 


0/ ( 2.0-FF)  I 


45 


4’  CONTINUE 
IN  C AS  F OF  ATH.GT.AO 

IF ( T.OT. TR)  GO  TO  49 
A 1 I ) =A  0 A''  ( ( 2. 0-.E ::  9 ) /2. 


_Q1 


A ( I ) = A ( I )«•(  ( 2. 0-E  FD)  / 2 .OKCDSW  E F D7Tr WWWWIUT- 
l M PI /2  .OM  * < EF3/2  .0) 

IF  ( A ( I ).LT.  0.00001)  33  TO  45 
_A  1X1=  A ( U1XI2.0/T  2.0-  E.E.D  )J 


C 

GO  TO  48 
C 

4 9 A (I  ) -A  THA  1 { ( . 0-E  F Ml  .0) 

A(  I )=A(  I)*{  ( 2.0-EF  )/?.  0)  *:*Q;  M FOREMANS'  M (BE-EF  ) /BE  ) * ( T-TB) 
1 * (PI /?  .OH*  ( EF/2.0) 

IF(  A(  I ).  LT.  0.00031  ) 30  TO  45 

A ( 1 1-  A ( 1 1 5 1 2 .0/1  2,0  - E F LI _____ 

00  TO  43 
4 5 A ( I ) =10000.  0 

0 

4.8 IF  ( FSF.E0.3. 01 AS  =1033  .3  

I F ( ABS(AU)  ) .GT.AS)  30  Tn  15 
8 CONTINUE 
GO  TO  46 

1-5 KK  = I _____ 

00  60  I=KK,NPT 
6 3 A ( I)  =A  S 
46  CONTINUE 

C »*£ All  ).  CAL.  END  Mi 

I F ( ATH.GT.AC ) GO  TO  553 

TC=2.0/(  2.0-FF  K ( AC'  '( { 2. 0-=f  ) /2 . 0 ) -ATI  * M ( 2 . 0-  E F 1/2.0)) 

TC  = TC/C/Q*MEF  /BE  )/ANS*M  ( BF-FF  ) /3  E ) / ( 3 T / 2 . 0)  * * ( FF  /2.0) 
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TC  = TC  «-TB 
GO  TO  364 

36  3 TC  = 2. 0/(  2.0-EFDd  ( AC**  (I2.D-FFD1/2  .0  1-  A3**  ( (2  .O-EFDI/2  .0  M 
TC=TC/CDS/3  **(  FF3/3r  )/ANSf'5'  ( I 3E-E=D)  / BE)  / ( * I /2. 0)  KM  EFO/2 .0) 
-36  4 WRITE(6,800)  S Ml  1 .CMP  jJJ  , S B ( 2 1 . CM P ( 2 1 , S R ( 3 ) # CMP  ( 3 ) , ZBB  , 3 . AN  S 


800  "FORMAT  ( !h  ,7  , • S B 1 , S3  2 ",  SB3  , 7 BB , 0,  ANS  * , / , 

131  El  2. 4,  • ( ' , F4.1,  ' ) *1  » 3F12.41 
WRITE  (6,801  1 ATH,T3,AC  ,TC,CDS,  FS  F 
_&0_L.£1RMXUH_^7,J A TH  , T 3 , A ; , TC  P S.  = » . / , 6 El  2 . 4 1 


C 

C RESIDUAL  STREM3TH  RATE  I FSF=1.0(FAIL  SAFE  I 
DO  18  1=1,  MPT 

IFlFSF.FQ.l  .0)  03  TO  189  


=0  • 0 ( NO  FAIL  SAFEI  1 


I F ( All  1.LE.AC1  00  TO  183 
R (n  = FLOAT{  KC  )»S3RT<2. 0/P  I /A  (I  ) 1/R  0 
00  TO  18 

183  R ( T 1 = 1 .0 _ 


GO  TO  18 
139  SW  = A8S  ((  A(I)-AO)/  ( AS-AOM 

RU  1 = 1 .0-C  l.0-ZETA»*S3PT(  $W) 

! 3 rONTTNIlF 


r 

c 


$$*  CAL  . OF  H $$$ 
nD  115  M= 1 , MPT 
.HMA  =Q_^Q_ 


OMN=R{  N) 
no  116  1=1,3 

IF  ( CMP  { I 1.  EG. 0.0)  sgc<  ii=sgs  m 
_VXGIL=SG_C(  I I /XO 


AM  FT  A ( I)  =GMM*AMUO/SG;(  I 1-XO/SGCU  1 
RR  ( I ) = V0*GMN*-AMJ3/SGCl  I 1 
ARS0  = ANETA(  I 1/ SORT  (2.0  1/RRU) 

.A  R.3  1 = - AM  FT  A ( LL±1U  I)*»2/?.Q 


A R32  = ( ANETA  ( I I — R R { I )M  2)  /SORT  (2.0  1 /R"RTFT 
ARS3=ANETA(  I)/RRin/S3RT<  2.  01 
ARG4=1  . 0 /SORT!  1.0  l-RR(  I 1**  21 
ARG5=-0»  5*  f 1 .-0/V X ( I 1 * » 2*  1 . 0 / VO* » 2 1 


ARG  6=0 .5  * l RR  (I)  / 7 X C I ) * 1.0/VO  1;;^2/(  1. 0+RR(  I ) ¥¥TS 
ARG  7 =ANET  A ( I 1 / SQR T ( 2 . 0 1 /RR ( I) / SORT ( 1 .0 *R R{  I 1**2 1 
AR  38=1.0/ S3  RT(  2.0  1 /SO  R T(  R R { I 1**2 f l .01 

1 ■(  (FL-l.OI’-RRddVXU  ) fFL/RRI  T 1/ VX  (I  1- 1 . O/VO  1 

IFfCMPd  1.E0.1.0)  CALLHCOMt  H(  Ml  , A RG3  ,7R  G r,AR(T27ARG37A!rG4,_ 
1ARG5,ARG6,ARG7,ARG8,7P III  , AM ( I 1 1 
IF (CMP (I 1.EQ.0.0)  CALL  HS I N ( H (N 1 , ARGO, AR G1 , ARG2, A R 33 , ARG4, 

LA775-tA  R3.6-.A_  R.G  7 .A  3 3 8.  P PA  I ) , A N ( I ) . U P C (I)  1 __ 

HW  A = HW  A + H ( M 1 
115  CONTINUE 
HI  M 1 =H WA 

IF  (HIM  ) .GT.l  .0  ) H-LNJleJL  .0 


11  5 COMTINUE 
C nw  CAL  . OF  H FMD  ### 

CAUL  OSF (DT ,H,TH,NPT 1 
A 1 PH  1=  A1  PHA-  1 


DO  208  T = 1,MPT  “ 

T=FLOAT(  1-1  )T[)T 
TBET  A=  T/BFT  A 

- 20  3 W ( T 1 =AI  PHA/8ET  AT  BET  A»  :A.  PHI  E X ? ( - T3  FT  A * » AL  PHA  1 

WR I TE(  6,10101 

1010  FORMAT ( 1 H , / , * I , TI ME , A , R , H, TH , W * , / 1 
C**MF  MEED  MORF  ACCURATE  PJTVOO  80  1 = 1 , NPT , B IGGER  •*'** 


-61- 


no  BO  I = 1 .NPT.RO 
T = R L 0 A T ( I - 1 ) •-  0 T 

WRTTE(6f  1023  ) T t T , A ( I)  , R ( I ) , H ( I ) , T H(  I)  , 4 ( I ) 

1020  FORMAT  { 1 H , I A , F 8 . 1 , ? F 1 1 . '+  ,3  F 1 0 . 3 ) 

3 0 £0  M TIM  UE ...  . . 

WR!TF(6,600)  NN 
600  FORMATUH  , *NN='  , T 5 I 
no  77  1=1, MPT 

01.1  MY  ( I ) -H(  T ) ... __ 

7 7 Z 3 M T I M Ur 
H0=OUMY(n 

c. 



C •!  b n v " «t!  i:  :A’  it-  if  ! 3*  S'  f\|  Q T \J  S P R 0 T I 0 M CASE  *•  >„■  V v }"  * t *■'  f-  t * < *’  ? *■  * F 

NN  = N 0 T - 1 

TI  =TT/FLOAT(  NNf  l) 

U.Q.T.  =_I  JZD..T. 

NNN=NN*1 
no  in  1 = 1, MNN 
TTO=TI AFLOAT < I ) 

IT.  =N  DJJOil 

DO  ?0  K=  1,1  T 
T=  F L 0 AT  ( K — 1 P DT 
\H  = I Tfl-K 

X.=  - Jjt  ±L0-  ,IHi  M HJ 

FX=EFXP(X) 

niJMY(  K )=W(  K ) fW(K)  fr  EX 
20  COMTIMUF 

-04  LJ O.SF.XOI.^JOUJ4Y.^-!1UMY-^J.ILL 

X = - TT Or“H  0-(  TTO/Dr  TAP  * ALPH  A 
XX=EFXP(X) 

P(  I > =-XX-D!)MY(  IT  ) 

I F ( P.  { I J.»L  T • 1 -0E-201  ?.(  1 ) =1.  0"  -20  

WRT  TE(  6,683  II  , IT  ,TTO,X  ,XX  ,OUMY<TT  ) , P(  I) 

6R3  FORMATdH  , I 3 , 1 6 , F 1 0.  1 ,4C  1 5 . 7 ) 

10  COMTIMUF 

0D_._j-0 30.00 : 

* * • * * iff*  *.  <•  <*  *• » * * If  >r «.  I N S P E C T ! OM  f ASF  ;•**■>  *»•  FM  wM  ? A ?• 

100  COMTIMUF 
NNM  =MMf 1 

. T I = I T/FL  OAT  (MM  FI  I 

MOT  =T  I/DT 
I T =NDT  Fl 

00  2000  J=1 , MNN 

XI  a=  IJJiELO  A-T-IJ4 

r * « p s ( J > t-  * 

DO  200  K = 1 , 1 T 
T = DT^F  LOAT  ( K -1  ) 

NW  = ( J-  IMMDTfK 

X=-(T  4-FLOAT  ( J-l)  ? T I UHO-TM!  MH  ) 

F X=  EE  XP( X) 

hum y (.<  j=w(  mw  )fw(  Ma*  hex _ 

200  COMTIMUF 

CALL  QSF  (DT,DUMY,DUMYf  IT» 

X = — C FL  0AT(J-1P  TT/RFTA  )-  ^ ALPHA 

Y=-1FL0AT(  JILT  1/3  ETA  P;;L  AL  PHA-Fl  .0  AT  (j  PH  I/  .HO 

EX  = EEXP(  X) 

FY  = E FX  P{  Y ) 

PS(  J) =EX-EV-OUMY(  IT) 


IF(J.NE.l)  GO  TO  210 
P(  J>  =PS(  Jf 
GO  TO  220 
21  0 CONTINUE 

r fc  it  k »;  t.  fc  %■*  11  ± 01  . J fr  j*  ft  »i  * t-  t *: 


WP  T T E ( 6,  1009  ) 

1009  FORM  AT  (1H  , 3X  » 1 QT  J ( I)  = * ) 
JM  = J-l 

-D-O  1000  4 = 1 tJ-M 


DO  250  K=1  * I T 
T=0  T^F  LO  A T(  K-l  ) 

AM={  1-1  I^NOT  +'< 

-22  a TuguagAT  u -ixt  u-u_l 


C * *•  «- 


T-  R M 1 
NA  = IT+ 1-K 
NH=N  A 
.MJU24- 


X=2Z-TH(  NH ) -RK ( N<  ) 
EX=EEXP(X> 

A X = A { N A I 

-XER  M .1=  — F-<  AX  »-  A1  ,A?,P.M,  J LTEX 


C TERM  2 

FSM=  l. 0 
JII=J-I 


NA=KK*  NOT*-  2-K 
AX  = A ( N A) 

FSM=FSMr  (1.0-F(AXf  Ai,A2,RM,un 
-MH=  ( J - I>-Hr.NnT>2-K 


350 


b * 


X=  Z Z-T  Ht  NH> 
EX=  EFX  P( X ) 
TERM2=-F  SM*E  X 

T-ER-M.-.-3-XjX- 


TER  M3  =0. 0 
I I = J-I 

I F(  I I.LT  .2  ) GO  TO  480 
JUL-AjOO  XK-2-tuXI - 


FSM=l.  0 
KKM  = K<-1 
00  460  M=1  » KKM 
NA  = Mi  MOI>  2-K 


-AAO 


AX  = A ( N A) 

X-F ( AX  »A1 , A2 ,RM,U> 

FS  = 1.0-F{AX, A1,A2 ,RM,U ) 
PCM=FRM»P<; 


NA  = KK*  NO  T*  2 -K 
AX=A(NA J 
MH  = N A 

±!K = J - 1 -KKJ-JL 


X = Z Z-T  H(  NH I -RK  C N < > 

FX=EEX P( X ) 

43  0 T=RM3  = TERM3-FSMSP  { AX , A l , A?  » R M , U I * E X 

_ 43-0 1XJMY(K  )=  CTFRMU:Tr  RM2XI.FR V13  HiaU.SLkU 

250  CONTINUE 

CALL  Q SF { DT  » DUMY, D UM Y» IT) 

QT  J ( I I=OUwT  ( IT  ) 

■IPOD  CON-TJ-NDE 


TOT  J = 0 .0 
DO  490  1=1, JM 
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mmunuu 


I = ( P(  J).  LT.  1.  OF.-?  01  P(  J>  =1  . 0F-20 
2n  PK(J|  = P(JI  + P(J)?,2/2.'l+MJ)v'3/3.0 

2000  CDNTI MUE 

C 

300  3 rn\!T  INIJF 

WR  I T E ( 6 » 59  3 ) 

59  8 FORMAI-tlM.  , /.,± I ) , T P(  I ) * . /) 

00  321  1=1,  NNM 

TTO=TI*FLOAT(  I ) 

IF(P(  I l.GT.0.0001  » S3  TO  324 

T P ( T)=FI.  OAT  ( NFMP(  !)  -FLOAT  jjj£  )* 

SO  TO  325 

324  IF(t>(  D.GT.0.9)  GO  TO  326 

tp(  n = i.o-(  i .o-p(  n mf 

GO  _J0  375  __ 

32  6 TP  C I»=  1.0 

325  CONTINUE 

WRTTEl 6,599)  I ,TTO,P(I ) ,TP(II 


321  CONTINUE 
GO  TO  81 
32?  STOP 
_END 

SUBROUTINE  TaTMTT , AROT AP"I7AP2 , AR3  , AR4  , ARE  , 
1 AR6  , AR7,  ARB,  P,  A\| ) 

C CAL.  FOR  COMPOSITE  GAUSSIAN 

M- 0.5*AN*PM  M .0-  E RF(  A ROD  *E  XP  ( AR1  H ( 1.  3 «•  E 

RETURN 

END 

SUBROUTINE  H SI  NM  , ARO,  API  , AR2  , AR3  , AR4  , AR  5 , 

LAB. 6 , A_RIj. ARB_lP,  AN i UPCO 

C CAL . FOR  SINGLE  GAUSSI AH 

IF  ( UPC  .EO.  1,0)  SO  TO  10 
C COMPLETE  SINGLE  GAUSSIAN 

H = 0.  5»AN»PM  (1  .O-ERM  AR3)  ) »AR4 

1 E X 0 { A R 5 4-A  3 6 M'  ( 1 . 0 t-F  R : ( A R 7 ) ) ) 


OATIN  F-  1 > L P(  I )«»2/2.0 


1.  3 » ERF  ( AR21  ) 


(ART*  )) 


GO  TO  25 

U^ER  CUT  SINGLE  GAUSSIAN 

_L3  H=  0 . 5f  AN  X’P  i ( (1.0-FRF(AR3)  ) <-  AR  4 

1 < : E X P ( A R5> A P 6 > > { E R F ( A R 7 M-ERF  ( A R 8 VI  i 
25  RETURN 
FND 

_EiMC  T1.DN  EEX  ELARSl 

A=0. 01 
R= -5  0 . 0 
C=1  . 0E-2  0 

IEJ.A B SOAR G!  . LT,:  )_S0  JO  .5  

I FI  ABS  (ARgI.LT  .A)  GO  TO  10 
I F ( ARC,  .L  T.B  ) SO  TO  30 
EEXP=EXP  ( ARG  J-1.0 

_GJ1  JLQ 20 

5 E E XP  = A RS 
GO  TO  20 

10  EEXP  = ARO*ARG-'-ARG/  2 .O  + ARG*  ARG’  A RG/  6.0 

_ GO  10  _2D_ 

30  EFX  P= - 1 . 6 
20  CONTINUE 
RETURN 
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FUNCTION  F(X,Al,A2tRM,U» 

I F ( X.LT.A1)  GO  TO  10 
TF(X.GT.A2>  GO  TO  20 
_E=Xt  X-Al  I /I  A 2 - AUJJLL3.N 


FUMC  TION  F S ( X) 

FS=l .0  -F (X) 

RETURN 

-EMO 

FUNCTION  ER  F { X ) 
IF(ABS(X).OT .10.0  ) 
X=X*  1.  414214 


GO  TO  3 


0=0.398942  * EX  P( -X*X/ l .0  ) 

ER==2.^(  0.  5 -Of  T-M  { ((  1.  332  74T  T-l.  821255)  * TH  .731473)  * T 
1-0.356563  ) ^T+0.3 19331  M 

IFTXl  1^2+2 

ERF=-ERP 
GO  TO  2 
■ ERF=1.C 

T.F  fX.  IT. O.OI  FRPs-FR: sag 

RETURN 

END 

SUBROUTINE  0 SF ( H , Y , 7 , NDI M ) 


di nens  ion  r<ii«zm 

HT=  . ■2-5'a^'a'a'a*  U 


I F ( NDIM-5)7»  8,  1 

NO  I N IS  GREATER  THAN  5.  PREPARATIONS  OF  INTEGRATION  LOOP 

— SXMX=Y  (2  H-Y  (2) 

SUN  1 = SUM  1 *•  SUHi 

SUN1  = HT*(Y(  1 H-SUNll-YI  3 l> 

AUX1  =Y  (4  )*Y  ( 4) 

AJX  1 = SUM  U-HT*  ( Y(3)  +-A  JX14-YI  5 M 

AUX2  =H  T*  (Y(l  H-3.3  7 5MY  (?)  «-Y  ! 5 1 » *2.625*<f  ( 3 > +Y(  4)  )«■/<  6)1 
SIJN  2=Y  ( 5 ) <■  Y ( 5) 

_ S-IJM2=-SIIM  ZfeSll 
SUN2=AUX2-HT^(  Y{  ^ H-SUN24-X  {b  )) 

Z(1)=0. 

AJX  = Y(3M-Y(3  > 

A! I X = A I 


Z(  2 ) = SUM2-HT*{  Y(?) +AUX +Y(  4)  ) 
Z(3  ) =SUM1 
ZI  4 ) = SUM2 

-A! 


INTEGRATION  LOOP 
DO  4 1=7 » NDI M, ? 


SUM  1 = A UX 
SJM2=AUX; 

A'J  X 1 = Y ( I -1  > <-Y(  I-l  J 
AUX1  = AUX  1+AUX1 

JtUX  1 = S UM  1 f H T ft  ( Y ( T -2  1 f A LLX 1 +Y  ( T ) ) 


7(1-?)  =SUM1 
I F ( I-NDIM  )3»  6,  6 
3 AUX2  =Y  ( I I Y C I » 


AUX2=SUM2+HT*  ( Y(  I - l ) * A UX  2 +Y  ( I f 1)  ) 

4 Z ( I - 1 ) =SUM? 

5 Z(M0IM-1)=AUX1 


RETURN 
7 ( M D I M - 1 ) = SUM2 
Z(  MDTM  ) = AUX  1 
-RE  T-URN 


C 

c 


END  OF  INTE  GR  A T I D N LOOP 
7 IF(NDTM-3)12,11,8 


NDI M IS  EQJAL  TO  \ OR  5 

SUM  2 = l .125*  H T <'  ( Y ( U*Y(  2H- Y(2)  «■  Y(  21  *Y(  3M- Y ( 3M-Y(  3 U- Y<  * M 
SJM1=Y(2 )+Y (2) 

SUM1=SUMI»SUM1 


SUMl  = HTMY(  I )+SUMlfY{  3 M 
Z(l)  =0. 

AIJX1  = Y(3)*Y(  3) 


7.(2)  =SUM2-HTt(Y(?  ) fAJXlfY  ('+  ) ) 
IF  C NDIM-51 10,9,9 
A JX 1 = Y (4)fY(4) 

AUX  JL=  AUX  m u x 1 


Z( 5 ) = SUM1+HT* (Y (3  H-AJX  H-Y( 5 1) 
10  Z ( 3 ) = S UM  1 
Z(  4)  = $UM2 

RETURN - 


C 

C 


L 2 

//CO. 

-123. 

/* 


NDIM  IS  EOJAL  TO  3 

11  SUM1  =HT*  (1.2  5*Y(1)  *Y(2  M-Y  (2  )-.25*Y  (3  )> 

SUM  ?=Y  ( 7 1 f.Y  ( 2J 

SUM2=SUM2+SUM2 
Z ( 3 ) =HT*  ( Y(  1 ) *■  SUM2  f Y ( 3 M 
Z(  1 ) =0. 

Z ( 2 )-SUM 1 


RETURN 
END 

SYSIN  DD  *■ 
0 
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APPENDIX  C 


«f5(F.sS§Wh&&  *,ncs= 


3LKSI l E=  3440 


ABE  L2 

SENSITIVITY  STUDY  IS  AVAILABLE 
OE L T-A  K T H RESHOLD  IS-AVAI.  ABL  E 


HREF  LOADS  (COMPOSIT  OR  SINGLE  GAUSSIAN)  ARE  AVAILABLE 
ANY  fOTAL  FLIGHT  HOUR  IS  AVAILABLE 
DIMENSION  AA(10)  , AAAI  10  ),  3AB(  10  ) 

f)I  m^sION  - A-A  AAM-t-04 

DIMFNS  ION  A ( 1501 ) , R ( 1 5 01 ) , -(( 1501 > , TH ( 1501 ) 

DIMENSION  PP(3),SGC(3),RR(3),  AN  ETA  (3) 

DIMENSION  CMP(  3)  ,SGS(3),VX(31  *AN(3  ) , UPC  < 3 ) , S B(  3 ) 

— OIMFNS  ION  0 ( 126) , GAf 1?  6,  15  ) , G TH(  1 2-6rl5)  , I T0-LP5) 

DIMENSION  G ME ( 1 2 S ) 

COMMON  Al,A?»RM,U,NT,NPT,DA,NF  ,PME 
COMMOM/QWE/  ALPHA,  BELA  , EFD,  CDS  , 0 , A NS  , BE 
'COMMON  /POI-/-e-A-l-rMC-TNR — — — 


a?  &'  ii  & $ & It  & ^ $ # W 1*  1*  ¥ & 

is  fe.atJ^Lle— 


31  CONTINUE 

p4s*i  *£*:**.'  $ ft:  <;**!*♦•  if  L .tW  A ♦sf'F  *****#i>:** 

C * * **  **********  MEAN  VALUES  )f  **«-**  ' ‘ * 

FRACTURE — DATA — 

FF  =2 . 5 
8E=4.0 
DK  T H= 1 .5 

- — EED-KL.O — 

C = 3 . OE -7 
KC  = 6 0 

F $ $ *i  u dr  * * >>  * * * -it  * 1 * * * INITIAL  CRACK  DATA  * # *»?  t *■  * * * **  *■  tt  f if  * * 
„1 5 

A0=0 .04 
AU  =A  0 
AL  = 0.0  15 


A(  PHA=4.0 


RE  T A =3  0000. 0 

f **  *.**#•*  ******  I*.  * LOAD  DATA  ***  * *«?•  k%h  * **  ?»  <t?  *.  $ *{?  *■ F k 

C CMP  ( = 1.0  COMPOSIT  GAJSSIAN  ,=0.0  SINGLE  ONE) 

CMP(  1)  =1 . 0 


c.  uJ: 


C EL 


CM  P ( 2 ) = 1 .0 
CMP  (3  ) =1  .0 

(=1.0  UPPER  CUT  SINGLE  GAUSSIAN 

UPC(  1 > = 0-.0 

UPC ( 2)  =0  .0 

ijp:(  3)  =0.0 

UPPER  CUT  LEVEL  OF  SINGLE  GAUSSIAN  (EL^XO) 
E-tr-S-s-e- 

G=  10. 0 
AA(4)= 115.0 
AA A ( 4)  =115.  0 

A A A Af  44®  00 . 0 — 

BAB (4) =3 .0 

C AN  ( I ) NUMBER  OF  LOAD  CYCLES  PER  HOUR 
AN(  1 )=  600.0 


, =0.0  COMPLEr  E ONE) 


AN  f-2  1=600.0- 

A N*  ( 3 ) = 60 . 0 
ANA  = 0 . 5 
PP( 1 ) =0. 4975 

ppf  21  = 0.  002-5 — 

PP  ( 3 ) = 0.  5 

C S3 C STANDARD  DEVIATION  OF  COMPOSITE  GAUSSIAN 
SGC(1)=0 .07 *G 
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m 

f S3S  STANDARD  DEVIATION  OF  SINGLE  GAUSSIAN 
S3  S I l ) = 1 . 5*  G 

SGSL  2 ) -1. 5iG 

$G  S ( 3 ) =1  .5-'  G 
XO  = 1 . 0'?'G 
ZR=1  .5 

f At*  *■*  -fr********-  **•  MAT  - RT-At--  DA  T-A--*±*^  *•**•*•*> 

V0=  0 .0  56 
AMU0=5.7*G 
R 0=  5. 7^G 

. £ Jt  |fc  * *.4  ^4. * Ifc^ - E-O-A  TA «:i:w  ft  * i •■•  ft  fe  r <c  * A ft  A Jtfr  A k t )t — 

FSF  =1 . 0 
7 E TA  =0 . 43 
AS  = 7 .0 


■A*-! 


iufct  Jt: 


-Vi^  -4-  * - 


r 4 4j»  - 4 i:  At  £ * =V  Ask  4 A ~ ■(<  * A F A 

C TT  TOTAL  FLIGHT  HOUR 
TT= l 50  00 . 0 

C * * ► i<  * A * k A k if  ■*■  $ * * i * INSPECTION  DATA  *■<  ***:}•*•**  ***#«#**#** 

- A 1 = 0.0  2 

A2  = 0 . 3 
PM=0. 1 25 
U=  1 . 0 

C NN  NOMRER  OF  TNSPEG-r  IONJ5 


C TT/(NNH»  MUST  RE  INTEGER 
t (TT  / CNN*-  II  )/(TT/(NPT-l))  MUST  BF  INTEGER 
NN=5 

r :-■■■  I — — A OF  I N SP  EG-  .— GT  . -5  — G A { — , 6 1 , GT-H  { r L ! .1 T J ^ { G 1 — MUS  T — BE — CHANGE  D 

C # OF  INSPPC.=NN  G A ( ,NNH),.3TH(  , NN  H ) , I TO  ( NN*- 1 ) 

(■+%*«  * :!.}C4*  \ >•:*-*  i * PROOF  TEST  DATA  *-■  cr*  vs--**  tfT-.’s  if  F* 

a;t=9rp.  o 

r AC  T - CUT-  LEV  EL  OF  -TNT  .CRACK  D ISIR  I3JT-I-0N 

C AGT.GT.AU  MEANS  NO  PROOF  TEST  PERFORMED 

CURSORY  I NS  RFC.  DATA  A:-****  i?:  MUH  4 IU 

C SR  Y = 0.0  -- - — — — 

NR  = 1 
NC  =4 

C A 1 = 0.  6 , , . 

-Q4.i-.Sc.  AwafcJ:  -W*..  ..  ..  s.  F = i ■ - ..  ll -I*:.-..  s.4L4ii.  - jtajL 

PI  =3.141  59 

C N’T  3ASIC  DATA  POINTS 

C EACH  DIMENSION  SIZE  0 F( A, R , H, TH, DUMY, A >MUST  3E  G.  = . NPT 

r (TT  / ( NP  T-l  -H  - MUST  RF  I NT  EG  ER 

MPT=  1501 

G NOT  DATA  POINTS  UNDER  NO  INSPECTION  CASE 

C FA  OH  DIMENSION  SIZF  OF(TP  P , d S , R K , C I J , CT  J ) MJ  S T BE  G.E.  NOI  AND  ( NN  *-1  ) 

C (T  T/N3I  ) - MUST-  BE  INTEGER  

r (T T/MOI) /(TT/ (NPT-1) ) MUST  RE  INTEGER 
NOI =30  . 

C ;>i>  }■  <■  ;*f:v4  **  -y  *■■■*- 1 * FLEET  SIZE  *-  * - -*■■■ ' -+c' x'  > ■ ■ 

- - - NF  =5  0 - 

C>?  Cf  •*•  )■> ■■■'■*"■■*•  ■ :H'5« f ' f *•?.  - "•  » • ' •'•  ' L » ’>  - f : 

Q -X-  nit  it  5:  ’/■  S'  *?■  *;■■*  A%V  ' v > • si:  * A v?.  4f  -*  z'  *•  - \ Y £ r u *•.  ^ . 

f*  U ^ ».  >:  »'■  ^ ;V  *?r  a> k+  . *:■  r-V‘J* o':  'c  *:  A- .*  ^ V r f-  ••  ,j.  . '■  !j'  : 

C SJB  PROGRAM  FOR  SENSITIVITY  STUDV 

RF  AD  (5,501)  XY»  YZ  1 » YZ  2,  Y7  3 
501  FORMAT (4F10. 0) 

WR I T E ( 6,  650  ) X Y,  YZ  1 , YZ  2 , YZ 3 


>■  T * f.  *•  V ts  * *•  i-  V i.  t (5  l: 


t r'  '4 f.:i'  * ; 

i'  jt  H ?>c : 


: 5 r i - 
:< ' j:  fc 

kit 


-68- 


650  FORMAT  (1H  ,///.*  XY  ,YZ  1 »Y  72  ,Y  Z3 » 
IF  ( XY.LE  .0.6  J GO  TO  32? 


I F ( X Y . F3  .1.0 > FP  = YZ 1 
IF  { XY  . EQ  .2 . 0 ) DKTH=Y  Zl 
I F(  XY.EQ.3.  n)  EF0=Y71 


IF(  XY.EQ  .4.0  ) C=YZ1 
IF(XY.EQ.5.0I  AO  = Y ZI 
I F < XY.E3 .5.  0>  AU=YZ1 
4F4XY  -.€0^6.04  -ANC4-1^Y-Z4- 


IF ( XY. E0.6. 0»  AN( 2 ) =YZ1 
IF  ( XY.EO  .7.0)  ANA  = YZ  I 
IFUY.EQ.8.0  ) AN { 3 ) =Y Z 1 
-4F44CV,-F0^9.-Q-1~-  VO~-Y  Z r 


IF(XY.EQ.IO.O) 
IF ( XY. EQ . 11 .0) 
IF{ XY.EQ. 12. 0) 
I-M  X Y-.  NE  .4.3  . 01 


Z E T A = Y Z 1 
A5  =Y  7.1 
KC.  =YZ  1 
00  TO — 7 24- 


PP  ( 1 1 =YZ1 
PP( 2)  = YZ2 
PP( 3 )=YZ3 

-424---OON4I-N-OF- 


I F ( XY.NE .14 .01 
SGCC 1 ) =YZ1 
SOM 2) =YZ2 
SO 04  34  —YZ3 


GO  TO  722 


722  CONTINUE 

IFCXY.NE.15 .0) 
G = YZ  1 

— S-0C4-14— 0-.  0 7 - G - 


GO  TO  723 


SG  C ( 2 ) =0 . 1 R * G 
SGGl 3) =0. 1* G 
X0=1 ,0*G 
R0=5.7*G 


AMU0=5  .7*G 
723  CONTINUE 

IF  { XY ,NE .16 .0) 
-4MU0^734 


GO  TO  7 24 


R 0=  YZ  1 
7?  4 CONTINUE 

IF(XY. EQ.17 .0) 
4-F  ( X Y.  €-Q  .4  &-.-Q4 


IFtXY. EQ.19.0) 
IF(  XY.  EQ. 20.01 
I F ( XY.EQ  .21.01 
-4F40CY-.  EQ  . 2 1 .-04- 


AL  PHA  = Y Zl 
44=44=47-4— 


FS  F= YZ  1 
BE  =YZ1 
A\( 4)=YZ1 


IF  [ XY.FQ.22.  0) 
IF( XY.EQ .23 .01 
IF(XY.NE.24.0) 
SGS4  4)  =Y  Zl 


AAAA(M  =Y7  1 
ZR  = YZ  1 
GO  TO  725 


- SGS  ( 2 )=Y  Z2 
SGS  ( 3 > =YZ3 
72  5 CONTINUE 

IF ( X Y . NE .25 . 01 


GO:- TO  ? 26 


CMP ( 1 1 =YZ1 
CMP ( 2) =YZ2 
CM  P ( 3 ) =Y  Z3 
72  6 CONTINUE 


IF(  XY.EQ  .26 .0)  A1  = YZ1 
IF(XY.  EQ.27.0)  A2=YZ1 
I F < XY.  EQ.28.  0)  RM  = YZ1 


4F15.7I 
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rSo 


Iflxv: 


XY, 


■ :§3 


EQ 
EQ 

IF  C XY.  NE  .31. 
UPC  l 1)=YZ1 
HP 6 F 2)  =Y72 
UPC( 3) =YZ  3 
7>7  CONTINUE 

I F ( XY. EQ  .32  . 

- I PIXY. EQ  .33. 

IF(XY.EQ.34. 
C SJ  3 PROGRAM  FOR 

«•  sfe ;**■  *f;  :#:<  fr*': 

JKL  ®0 


8!  HsZfh 


0) 


H TO  727 


0)  EL  = Y71 

0)  437  = 77  1 

0)  C A l =Y  7 1 

SENSITIVITY  STUDY  ***END 


t*  i 


V t'-f 


f:  & 


■ k js  fc  jf  ■ F r>!  $ k tf.  ^ t’  k: : 4 sf:  & $ fc  jt  }e 


KB  E =BE 
I RE  = 8 E 

IF  ( CMP  { 1 ) .EQ 
TPI  CMPU  )*-0 
IF ( CMP ( 2 > .EQ 
IF ( CMP (2 > .FQ 
IF  ( CMP(3).EQ 
I F F C M P 13  H E <3 


SR(  1>=AA(  TBE  ): 
-S8-F7rF^AA44F€-)- 


AB  ( I BE ) *?-S  GC(  1)  *K3E 
^SG$(l  d4  < BE 


.1  .3) 

-r0-r©~) . _ . 

.1.0)  SB  ( 2 )=AAA(  I BE  )*B  AB C I3=»*SSC<  2 ) < 3 E 

.0.0)  SBI2  ) = AAAU8E)  < S GS(  2 )**<  3 E 

.0.0)  SB(3>  =AAAA(IBE)  *SGS( 3 ) F*  K3E, 

.1,0 ) S B { 3 ) = AA  AA ( I B E ) * 3 A R ( 1 3 5 ) " ;S G GI  -3  F*  F < 3 E •• 


ZBB=(  ZR-G)  **  KBE 

Q = A N ( 1 ) t PP(  1 H SR  ( 1 d A N ( 2 ) f P P ( 2 ) * SB  ( 2 ) <- A N ( 3 M"  PP  ( 3 d S B ( 3 ) 
UAN  AT7.BB 

- AC;  ~2  *7  PF* ( Kf.AROd*? : — 

DT=TT/FLOAT<  NPT-1  ) 

ANS  = PP(1  )*AN  (1  dPP  (2)  * AN  ( 2)HP  ( 3 d AN(  3) 

Av)S=CYCLES  PER  HOUR 


IF(EF.E0.2.0  ) EF=  2 .01 
I F ( EFD  .EQ.2 . 0)  EFD=2.01 

A NT  H = A LOG  1 3 < C ) +EF  * ALOG  13  I DKTHF- 
C 

CDST  = ANTH-EF0^AL3G10(DKTH) 

CDS  =10  .3d  COST 

t:.±  

Q -frt*  k ilk  k r r5;  rs  y f.sk  ■ ’Vrl-.v-f'  •fc-Vskr  k if*: 

IPICSPY.EO.l  .0  ) NN  = NRMNR  + 1)“NC 
TF(NN.EO.O)  GO  TO  24 

NM1  = NN+1  - 

DO  20  I = 1 » N N 1 

20  ITOU  ) = I'MNPT-1  )/(  NNf  1 )H 
GO  TO  22 

2 4 r.  ONT  I N IF 

no  23  1=1,4 

23  TTO  (I)  = Id  NPT-1  )/  6+1 
NN  1 = 6 

2-2-  COMT  INUE  - - - - 

NRN1  =NBNd 
NT  = NBN  1 

CALL  DENST  ( 0 , N T,  AU , AL  , DA  ) 

IF  F ACT.GT.AU4  -PME-=  1.4 

TP ( ACT.GT.AU ) GO  TO  830 
NT  C=  ( AU-ACT  ) / DA 
MTT  = MT-NTC 

- no  831  I =NTT,NT- 

0(1 ) =0.0 
831  CONTINUE 

CALL  QSFIDA ,0,GME , NT) 
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PME  = 1 . O/GMF ( NT ) 

83  9 CONTINUE 

k fr  is  «£:!•<•??  5. < •ftjfc&g-  k%  £ Jts  » { %<■$■ 

F * #1  ! ' f 4f  {'  »:  4.  4 * * *V-  y.  * * 3^..  ,j.  4 1 f f : j*  C-  jj,  ■?■ 


F-TL^9999 

KILE=9999  MEANS  'CRACK' 
DO  31  KI  L=i  , NBN1 


TAKEN  OFF 


-e 


AO=  AU- FL  OAT ( KIL-1  I ADA 


AH  = { D KTH/  { Q/ANSM  Ml  . O/BE)  l**2*<2  .O/PI) 

TB=  2. 0 /{  2.0-EFD  )/:0S/3*M  Ep  D /B  E ) k A NS**5 M ( EFD-BE  )/BE ) 
TB  = TB/  (PI/2. 0)**(  EFD/2  .01 

-^3-=T9*  (A  Tf|H  ( (2.  O-ETD)  /2.-CH— A0»»  { {-2.0-EFD  )/2-r044 - 


,FP25U2«9rEF>/2.3»*C*3*«!(FF/BE»*ANS*-M(BE-EFI/BE) 

IM  P I /2.0J** ( EF  /?. 0) 


lf-  (PI/2.0  )«{  EFD/2  .0) 

DO  8 1=1, NPT 
T = FL  OAT ( 1-1  >*DT 

— FFTMBMsHThrD-) — &B-99-V3 

IN  CASE  OF  AT4.LT.A0 

AC  1 ) = AO*M  ( 2 .O-EF  )/2 .0  ) 

A(I»  =A(  I I *FC  T 

I FlTM  Obl  . L T .0  vOTXLO  14 -44-40-44- 

A(  I ) =A  ( I )**  12.0/  ( 2.0-EF)  ) 

GO  TO  98 
93  CONTINUE 

LF4-4i-S-4r4ft  f — GO — 49—49 


C IN  CASE  OF  ATH.GT.AO 

A ( I I = A0**(  (2  .0-FFD)/2  .0) 

A ( I ) = A ( I I fFCQD*  T 
-191  Ad4.  LT  . 9 *903 0 19-30  40-45 


Ad  )=A(I  )*M  2.  0/(  2.0-EFD)  ) 
GO  TO  98 


99  A{  I)  = ATH**-(  ( 2.0-EF  >/2.0) 

AC  I ) =A(I  )+ECQ*  (T-T  B) 

IF!  Ad).  LT.  0.00001)  30  TO  95 

Ad  HAT!)**  12.0/1-2  M3-E  FH 

GO  TO  98 
9 5 A{ I ) = 10000. 0 

- 9PT  FT4S4.F9  iO-.-m-A-SM.409  4- 


IF(  ABS(A<  I)  ) .GT.AS)  30  TO  15 
IF  ( KILE.  NE.  9999)  GO  TO  915 
8 CONTINUE 

G9-T-9  46 

- L 5 KK  = I 

DO  60  I = KK,NPT 
50  A ( I t=AS 
- 96-C-ON-T4NATE- 


GO  TO  91  6 
915  DO  917  1=2,  NPT 
917  ACT)  =A  ( 1 ) 

- — 9L444NT40U9 

C un  Ad)  C AL  . END  ### 

IFCKILE.NE.9999)  GO  TO  930 
A33=A0*1.2 
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C 


IF(A(NPT).LT.ARC)  KILE  = KILH 
43  3 CONTINUE 

I F ( ATH.GT.AC ) GO  TO  3S3  ^ . 

t:=?.0/<  ?.0-FF  ) ( AC*M  ( 2.3-EF  > /2 .0) -ATHf  M (?  .0-EF  ) / 2 . 0 ) ) 

T C = T C / C/  6* M E F / R E ) / A M S i M ( a E- EE  4/0047(0  1/  2 -a4-i^4-E-E -/  2~  Q4 — 
TC  = TCf  TB 
GO  TO  36  4 

33  3 TC  = 2 .n/(  2 .4-EFO)  ' { ACC  * ( ( 2 .0 -EFD  1/  2 .0  I - AOv-M  (2.0-EF3)/2.0)l 

i:  = TC/CD  S / a-i  * ( FE  D /B  E4  /ANSi--  ( 4Sfe=*€4=  Ol/.ai-i  /4P- I /^a4-fc^4£4.Q42-*&4- 

36  4 CONTINUE 
833 
831 


CONTINUE 

CONTINUE 


( FSF=1.0(FAIL  SAFE) 


C 

c 


RE  S IDUAL  S TREMGTH  RAF  E 
03  431  I = 1 » N PT 
I P ( F SF  .EO.l  .0)  CO  T3  189 

Iff  A(  14.  LE.  AC)  GO  TO  183  

R(  I»=FLOAT<  KC)  $S3RT{  2.  O/PI  / A ( I ) ) /RO 
GO  TO  18 
R ( I ) = 1 . 0 

GO TO  18— 

SW=ABS  ( ( A ( I ) — A 0 ) / ( AS-AO)  ) 

R ( T ) = 1 . 0-  ( 1 . 0-  ZE  T A ) s S 3 R T ( SW ) 

IF(  K IL  E.NE  .9990)  r,0  T3  418 

40  CONTINUE  

GO  TO  419 
00  420  1 = 2, N PT 
R ( I ) =R  ( 1 ) 

-CONTINUE — 


= 0 .0(  N 3 FAIL  SAFE  ) 


188 
139 
1 3 


SGC  ( I ) = S3  S( I ) 


413 
4?  3 
41  9 

CAL.  OF  H 11$ 
no  115  N=l,NPT 
HWA  = 0 . 0 
GMN=  R ( N) 

00  116  I =1 ,3 
T E { CMP  (I  ).c0.0.0) 

—  VX  < I )=  SGC(  I 47X  0 — < 

A M E T A ( I ) =GMN-!'AMUP/SGC(  I ) -X3/SGCI  I ) 

RR  ( I ) = VO  4 GM  N A M J 1 / S 3 1 ( I ) 

I F ( C M P ( I ).NE.l  .3)  GO  T 0 21 

ARGO  = ANFTAf  I )/SPR  T4  ?.0>/RRU4 

ARG  1 =-ANET A ( I ) <-Rp  ( I ) **  ?/?  .3 

AR32=(  ANFTA  ( I )-R.R  { I ) * F 2)  / SORT  (2 .0  ) /RR  ( I ) 

IF  ( CMP  (I  ) „FQ  .1  .0  ) CAL.  HC3M  ( H(  N ) , AR  GO,  AR  S 1,  A R3  2 , A R 3 3 , AR3  4 , 

- 1 ARG  5,  A RG6  , A R G7  , AR  G-R  , P-P -(-14- , AM  (I  ) ) 

GO  TO  411 
21  CONTINUE 

AR3  3=ANFTA(  I ) / F R ( I )/SORT  (2.3  ) 

AR  3 4=  l . 0/S3  R T(  1 . 0 + RR  ( I ) - *24 — 

ARG5  = -0.5!l  ( 1 .0/ VX  ( I m?  + l .O/VO  ^<2  ) 

AR  G 6=  0 . 5C  ( P R (I  ) /VX  (I)  H.0/V0)^2/(  l . 0 »-RR  ( I ) ) 

AR  G7  = ANETA(  I ) / S3 1 T ( 2 . 0 ) / R R ( I ) / SQR  T ( 1 . Of  R R { I ) 721 

A R 38  =1 ,-0/SO R T ( 2 . 3 ) 7 S OFT  (ER  < I M ■ ? * 1 -04 

1 ( ( FL-  1.  0!  : RR(  I > / VX  (I)  4-FL/RR  ( I ) /W  (I)  -I  . 3/V3  ) 

IF  ( CMP  ( I ) . F 0 .0.0)  CA_.  HSIN(FHN),ARG0,AR31,AR32,AR33»AR34, 

1 ARG  5 , ARG  6 , AIR  G7  , AR  G8 , Po  (I  ) , AN  ( I ) , J P C ( I ) ) 

410  C ON  T I N IJF _ - — — — — — 

HWA=HWA*H(N  ) 

116  CONTINUE 
H(  4 ) =HWA 
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IF(H(N).GT.l .0)  HfN)=l .0 
IF(K  ILE.NE.9999)  GO  TO  421 
115  CONTINUE 
GO  TO  422 

4? 1 OO  423  1=2, NPT  

4?  3 H ( I ) =H  ( 1 ) 

42  2 CONTINUE 

*«  H CAL  . OF  H FMP  ### 

C ALL  OSFIOT ,H, TH » M PT f — 

KI P=NBNl-( < I L- 1) 

00  25  1=1, NN1 
G A { K I P ,1  ) = A ( ITO(I  ) } 
•GTH1KIP,  1 )=THMTDfl)> 


>5  CONT IMUF 
31  CONTINUE 

p <:|f  ?J-  :F:  jf-  Jy  V **:'  JS  !/?;**  7‘  5:  >;j?  £:  Uj'  V 

-e  subpou  time  c-in  • curs  a r y +ri  gorojs  —i-ns-pect  ions 


IFICSRY.EO.  1.0)  CALL  C INC  GA  ,GTH,0,  CAl  ,NC,  NR,  NN  , ITQ  > 
IF ( CSR  Y.  FQ  , 1 .0  ) GO  TO  81 

£ r*  S-  *t  $ ^ +#.*  ife*;  $:  fc  i-  *f>.  ** 

f SJ3P ou Tf NE  -)IM  PIG^OUS  1 NSPFCTION — 

CALL  H IN  (GA,  GTH,  0,  NN,  I TO  ) 

GO  TO  81 
322  STOP 
END 


SUB  ROUT  I NE  HCOM(H  , ART  , AR1  , AR?  , AR3  , AR4 » AR5  » 
1 AR6  » AR  7,  AR  8 , P ,AN  ) 

C CAL.  FOR  COMPOSITE  GA'JSSIAN 


-H  = 0.5*  AN*  P^-H  1.0-FRP(  APO>  )4-FXP{APl  )M 1 . 3 + ERF < A»? ))> 
RETURN 
FND 

SUBROUTINE  H SI  N(  H , ART  , AR1  , AR?  , AR3  , AR4,  AR5, 

1 AR 6 , AP  7,  ARB,  P,  AN,  UPC  } — — 

C CAL.  FOR  SINGLE  GAUSSIAN 

IF1UPC.E0.1. 0)  GO  TO  10 
C COMPLETE  SINGLE  GAUSSIAN 

H = 0.  5k  AN* Pi  < <1  .O-EPf  < A R3  » I *•  AR4- 

l'!'EXP  ( AR5f  A?  6)Y(  1 .0«-FR  = (AR7)  H 
GO  TO  25 

C UPPER  CUT  SINGLE  GAUSSIAN 

- 10  H=0 . 5PAN*P*(fl  . 0-  ERF(  AP3)-M-AR4 - - - - 

1TFXPI  AR5  4-AR6  )M  ERF{  AR7  M-E  RF  { AR8m 
25  RETURN 
FND 

FJNCT  I ON  ER  F (X) 


TF(  ABSIXi  .GT.10.0  ) GO  TO  3 
x=x*l. 414214 

T = 1 .0/ U .Of . 23 16 A 19  ABS  (X  ) > 

0=0.  39  894,2  * EXP{ -X1-  X/2 . 0)  - 

ERF  = 2.0(  0.5 -D*  T * ( ( ( ( 1.  33274*  T-  1.82  1256)  *T*1.  781478)  H 
1-0.356563  )*T+0.319381  )) 

I F ( X ) 1,2,2 

1 -ERF  = -€P.F - — 

GO  TO  2 
3 ER  = = 1.  0 

IF(X.LT.O.O)  ERF=-ERF 

2 RETURN — - - — - - - 

END 

SUBROUTINE  H IN  < A,  H , G,  N N,  I TO  ) 

01  ME  NS  I ON  A ( 126,1  5)  ,1(  125  ,15)  ,G(126)  , I T 0 U 2 5 ) , PX  ( 1 2 6 ) , P A ( 15  » 
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01  MENS  ION  TPM  15) 

COMMON)  A 1 , A 2 ,RM,U,NT»NPT  , 0 A ,NF  , P M E 
C N°  = 1 MEANS  SKI0  TO  WRITE 
N P = 1 

WRITE!  6,  600)  NN 
630  FORMAT  (1H  , • NN  =*  , 13  ) 

IFINN.NE.O)  vGO.TO  98T.,..  

S'  ♦;  i0c  *■.  * iji  * *.  ?r  4^  *•  -?■  !*•  s.  if.  f *•  v *•.  £ t ¥ * t ir 

f***  NO  INSPECTION  CASC  *^ 

N=6 

00-9-OOM--T-V-N 

00  901  J = 1 , N T 

931  PX!  J)  = OMEXP  (-H(J,  T)  ) * G ! J ) 

CALL  OSE  (DA  , PX  ,PX  ,MT) 

-  -----  P A ! H = PX  I NT  ) — — 

PA(  I)=PA(I)  TPME 

1 F { P A(  I)  .GT  .0.0031  ) 30  TQ.32V 

TPM  I ) =FLOA  TINE)  PAID  -FLOAT!  MF)  t F LOAT  ( NF  -1  KPA!  I ) M 2/2  .0 
no  TO  32  5 * — 

324  I F ( P A ( I)  .GT.0.9)  GO  TO  306 
TP  A ( I»=l  .0-  ( 1. 0-°  A ( I ) ) *TN!F 
GO  TO  325 

-  32  6 TP  A ! 14  =1.0  -— 

325  CONTINUE 
930  CONTINUE 

WRT  TE(  6,620) 

6°  3 FORMAT  UH  » RESULT  OF  N -0  I N S 0 F C . CASE*) — 

00  9 0?  I =1  ,N 

1 T=  ( ( N PT-  1)  / N ) If  l 

WRITE  (6,  621  ) IT,  PA!  I),  TP  A!  I ) 

-6? 1 "FORMAT  ! 1H  , 1 -I  T— 1 ,-T5t3-X-,  * PA  — * rEl-5 ■»-7-y-3-X-f-*-T-PA~  ’ ,~E-1§"»  74 — 

932  CONTINUE 
GO  TO  17 

C 

98  7 CONTINUE - - ----- 

F ft  if Jr  -jt  4 * !,):  I1!  •¥.  ► . :■  *s 1 ft  t * y * : ■"  -s  ;'” 

F -S'  fe  ( V A-  ij{  *!  4 It  if  £■;  ft  +;  & ft  ii,  •£  * f £ .»  p if-  s • •*  ft  >•'  r •>  v r ft  t ft 

F ft  & f INSPECTION  CASE  * 

IF!  NP  . EO  1 4 GO  - -T  O - 7 3-8  

WRITE!  6,699)  ( G( I I ,1=1  ,NT ) 

699  FORMAT! IH  , 5E1 5.7) 

DO  10  1=1, NT 

- - WRI  TE!  6,666  4 ! A ! I rM)  , M=!  ,6  4 -i - — - 

666  FORMAT (1H  , 6E15.7 ) 

13  CONTINUE 

00  2 0 1=  1,N  T 

-  - WR  IT  F!  6,  666  4 ! N!I,M)  , M = 1 , 6 ) 

?0  CONTINUE 
73  B CONTINUE 

C.  CALCULATE  THE  FIRST  PROB. 

00  -2 40  I = 1 , N T - - - - - - — 

2V  0 P X ( I ) = OM  FX°  ( -H  ( I , 1 ) M 3 ( I I 
CALL  OSE (DA, PX, P< , NT  ) 
p A ( 1 » =PXC  NT ) 

IfW.  EO.l)  GO  TO  88 B 
280  WRI  TE( 6,  170)  J ,PA ( 1) 

173  FORMAT!  /'PA‘,14,  * = *,F14.7) 
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83  R CONTINUE 

0#  iJ'  r S'  s-  *.  •>  4 :$* 

C CALCULATE  THE  ? - N PROB. 

01  100  J = 2,  N 

— no  140  - T=l , NT — 

C **  FIRST  TERM 
RH1  = H(  I,  1 ) 

RAi=A( I ,1) 

- RF1=FF  (RAIL 

TER1=0MEXP(-RH1) 

AR 1 =RF 1^  PA ( J -1  I 
T F R 2 = ( l.-TFR  l )«  AU 

T-F  R-M-l  =-TF  Rl-  MFP-? — ■ 

IF  (NP  .EO.  1)  SI  TO  2 60 

WRITE! 6, 270)  R HI , RA 1, AR 1, TER  1 , TER  2 ,TERM1 
270  FORM  AT  (6X  »E12.5,2F8.5,  2E 1 2 . 5,  El  4 • 7 ) 

**- T HTR  D — TFRM — — - 

251  FA  = 0 . 

JM1 = J-l 

DO  250  KK=1, JMl 

Ms  J — - • ■■  - 

R A3  = A ( I , M) 

RF3=  l.-FF{  R A 3) 

RH3  = H ( I,  M) 

— — — — 

TER3=RH4-RH3 

TFR4=RF3*((  1 .-OMEXP(-R  H3))*0MEXP{  -TER3)  fFA> 
IF  (NP  . EQ.i ) GO  TO  25  0 

35  0 FORM  AT"fr^X^2  El  2 ^ ^ 2FB  • 5^1 2 ^5^  E14  I 7^ ^ 

250  FA  = TER  4 
TERM3=  FA 

G--*T SECOND  FORM- — 

TE  RM  2=  0. 

IF(J.LE.2)G0  TO  200 
J2= J-2 

FA^TtO 

on  550  K = l,  J 2 
RH2  =H ( I,KH) 

R A 2 = A ( I,  KM) 

FA— AT  I tK4 — 

FA  = FA*  (l  .-FF  (RA)  ) 

RF2=FF ( R A 2 ) 

ARG2=1 .-GMEXP(-RH2 ) 

T F R M ~F  A-*-'FFF*  AFG?  h - P A ( J-K-14 

TERM2=TERM2  +TERM 
IF  (NP  . EO.  I ) GO  TO  55  0 

WRI TE ( 6,160)  RH7, RA2 ,R A, FA  ,RF2 ,ARG2,TERM 

IFF  FORMAT(2-X,FL-2T5r4TT.-5^E4F-^5-r^TT^-7T 

550  CONTINUE 

C SUM  OF  THE  ALL  TERMS 

211  TERM4=TERMl  4-TERM2  MFRM3 

p X ( I >-TKpMA  *G(  I ) 

IF  (NP  .EO.  1)  S3  TO  140 

WR ITE ( 6, 180 ) I , TERM4*  G ( I ) » 3 X(  I ) 

130  FORMAT (14, 1 * * * » 3 E 1 4 . 7 ) 

140  CONTINUE — 

CALL  QSF (DA, PX,PX,NT) 

PA ( J)  =PX ( NT ) 

IF(NP.EO.l)  GO  TO  100 
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m UHfUfc1™ 

no  ill  J=l  ,N 
P A ( J ) = PA  ( J ) * PM  F 

— I l l -FONT  INUF 

DO  887  J =1  » N 

IF  (f>  A(  J)  .3T.  0.  OOOl  I 3 0 TO  3 34 

TP  A ( J)  = FLOAT  (NF)  P A(J  ) -FLOAT  ( MF)  '-FLOAT  ( MF-1  ) ■«  A ( J)  S' f'  2/2.  0 
— ••• GO-  "TO— 33  5 

334  IF(PA( JI.GT.0.9)  GO  TO  336 
TP  A ( J)  =1  .0-(  1.0-PA  (J)  ) ;-"-MF 
GO  TO  335 

33  5-T3A+-Jl^l-.T> — — 

335  CONTINUE 
83  7 CONTINIJF 

1=  NT 

—  WRI  T frf  6*5101  — 

51  0 FORMAT  (4X,  1HT,  12*  , • A'  , 12  X ,»  H*  , 1 4 X , » PA  » , 11  X , • TP  A » ) 

WR  I T F { 6,  520  ) ( ITO  < J ) , A (t , J I , H(  I , J ) , P A ( J ) , TP A(  J ) , J = 1 , N ) 

520  FORMAT  (I  5 ,4E14.7> 

1-7 frO-N  T INUE — 

RFT'JRN 

END 

SUBROUTINE  C I N I A , H , G , C A1  , NC  , NR , NN , T TO  ) 

—  DIME  NS  I -ON  A 1126, 150,  HI  126,18-)  rG(  1 2-4-l-rTH^T-^^^--X14-244-^A14,8T 

DIMENSION  TP  A(  15) 

COMMON  A1*A2»RM,U>NT»NPT»0A  ,NF  ,PME 
C N°  = 1 MEANS  SKIP  TO  WRITE 

WRITE!  6 , 600  ) NN 
600  FORMAT  (1H  . *NN=*  , 13) 

IF(NP.FO.l)  GO  TO  733 

WR!TF(6,  699  > ( C-<  ! 1,1=  l ,NT  > 

699  FORMAT  ( 1 H ,5E15.7) 

DO  10  1=1,  NT 

WRITE(6,66S)  ( A ( I ,M  ) , M=  1 , 6 ) 

— 666  FORMATUH  6E  18.-7  ) 

10  CONTINUE 

DO  20  1=1, NT 

WRT  TE(  6,666)  ( H ( I , M ) ,M  = l t6) 

2-0 CONT  INUE 

73  8 CONTINUE 

C CALCULATE  THE  FIRST  PR03. 

DO  240  1=1, NT 

— 240  PX(  I+*-QM  EX  PT  -H  (1,11  1=0-11) 

CALL  OSF(OA,PX,PX,NT) 

PA  ( 1 ) = PX  (NT  ) 

J=1 

- IF(NP  .EQ.l  ) GO-  TO  888  — — 

230  WRT  TE ( 6 , 170 ) J , PA ( 1 ) 

170  F OR  MA  T ( /'PA‘,14,*  = ' ,F14.7i 
83  3 CONTINIJF 

A-isA-JL.  ■ — . — . — . — - — — 

C CALCULATE  THE  2 - N PR  03. 

00  LOO  J = 2, N 
00  140  1=1, NT 

C- *1 FIRST  TFPM  

RHl  =H  ( I,  1 ) 

P A 1 = A ( 1,1  ) 

NC 1=NC H 
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JD= JM i/NCl 
JOQ  = JO’tNCl 

IF{  JOO.FQ.  JM  1)  RF1=FF(RA1) 

I F{  .JQO  «MF  *JM  1 ) -c- l = F~.<  «A1  > 

T E R 1 =0  ME  X P ( - RH 1 ) 

A R 1 = R F It  PA  ( J-l) 

T E R 2 = ( 1.-TER1K  API 

— TE  RM1-TERH-T-ER? 

IF  (NP  .FQ.  U GO  TO  2 60 


27  0 
~e-+j« — 
26  0 


WRITE (6,  270)  RH1,  R Al, AR1 , T ER 1 , TER2 , TE RH1 
FOR  MAT(  6 X, El  2.  5 , > F 8 .5  , 2E1  2 . 5 , FI  4 . 7 ) 

-TH1RQ  -TERM  ~ - : — 


FA=0  . 

JM1  = J-  1 

DO  250  KK=1,JM1 

“M=J-Kk 

R A 3 = A { I , M) 

NC1  = NC  *1 
JMM  = J-M 

O Y—  J M M /NFrt 


JYY= JY^NCl 
I F { JYY  . EQ.  JMM)  RF  3 =1 . — FF  ( RA3  ) 
Ip{ JYY.NE.JMM)  RF  3 = 1.  - FC(  RA3  ) 


RH4=H  ( I » M+l.  ) 

TER  3=RH4-R4  3 

TERA=RF3=t((l  .-OMEXPl-R  H3  ) H'OMEXPI -TER  3 1 +F  A) 


350 

250 

C '*¥ 


WRITF(  6,  350  ) R H3,  R H4,  R A3 , R = 3 , TE R 3 , TF  R4 
FORMAT ( 6X,2E12.5 ,2F8. 5,F12.5, E14.7) 

F A = TER  4 

•TFR-M3*  FA  — , — 

SECOND  TERN 

T=RM2=  0.  . ■ ' 

IF!  J .LE.2  JGO  TO  200 
J2~d—2 


F A = 1 . 0 
DO  550  K=1,J2 
RH2=H{  I , KH  ) 

-R  A 2- ATO  fFv+  l ) — 
RA=  A ( I » K ) 
NC1=NCH 
JMK= J-K 

-d-T^d-MK-/NGi 


JTT= JT^NCl 
I F C JMK  .EQ.JTTI  F A = FA*{  1. -F  = { R A ) ) 
TF(  JMK  .NE.JTTJ  F A=F  At  ( l.-FCC  RA  )) 
dMKl»d- (OFHrF- — - — 


JU=JMK 1/NCI 
JUU  = JU^NCl 
Ic(  JUU.EQ. JMK1)  RF  2=FF  (RA2  ) 


AR3  2 =1  .-0MEXPI-RH2  ) 

T E R M =F  A!  R F 2 '•  AR G ? v PAC  J-K-l  ) 
TFRM2=TERM2+TERM 
+F-4-NP Ga^FF-3-5-0- 


££IIE.{£»160>  RH2»RA2»RA»FA,RF2»ARG2»TERM 
ISO  FORMAT  I2X,F12.5,4F8.5,  El 2. 5,  F 14. 7 \ 

550  CONTINUE 
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L 


c 


SUM  OF  THE  ALL  T ER  MS 
200  TrRM4=TERMl  f TERM?  HFRM3 
PX  ( I ) = TER  M4'G(  I ) 

IF  (NP  .EO.  1)  GO  TO  140 

WRI  TF(  6,  1801  I ,TERM4,G  { H , DX<  H — 

130  FORMAT (14, ' v^« ,3E14.7) 

IV 0 CONTINUE 

CALL  3SF  (DA.PX,DX,NT) 

PA  { J l = PX  ( NT  T - — - 

IF ( NP. EQ.l  > GO  TO  100 
WR  I TF  ( 6,  170 ) J , P A ( J ) 

IDO  CONTINUE 

0-0 11-1 J=l  r~M — — — — — 

P A { J ) = PA  ( J )>PMF 
111  CONTINUE 

DO  887  J =1  , N 

I F ( D A ( J 1 ,GT  . 0.  0 03  1 ) GO  T 0 -884 

TP A ( J ) =F LOAT (NFI ? P A ( J ) -FLOAT  (N F ) * FL 0 AT { N F- 1 J A { J 2/2.0 
GO  TO  335 

334  I F { PA( J ) .GT .0.0)  GO  TO  336 

— -TP  4 < JJ  = l .0-  ( 1 . 0-  ° A { J)  > 7?  fvip 

GO  TO  33  5 
33  6 TP  A ( J)  =1  .0 
33  5 CONTINUE 

83  7 CONT  INUF - - — - 

NR1 =NR H 

DO  26  1=1, NR1 

j=(Nc+mi 


l T&<  6 ,670  ) I,  PA  < J>  , I ,tp  A ( J > , !TP.(  J)  

670  FORMAT! 1 H ^ ( ' , I 2 , •)  =» , E 1 8 . 7 , 7 X , • ( • , I 2 , * ) = 


F15.7,l  OX,  15 ) 


>6  CONTINUE 
17  CONTINUE 

RETURN  - — — — - 

END 

FUNCTION  F c ( X ) 

COMMON  A 1 , A2 , RM,  J 

1 FIX.IE.A1I  Gn  TO  10  

I F ( X . G F . A 2 ) GO  TO  20 

FF  = U*  ( < X — A 1 ) / ( A2-  A l ) ) » *RM 

RETURN 

10  FF  = 0 .0  - - - 

RETURN 
20  FF  = 1 .0 
RETURN  ' 

SU3  ROUTINE  DEN  ST  ( G , N M, AU, AL , DA ) 

DIME  NS  TON  0(1) 

COMMON  /OWE/  ALPHA  , BET  A ,EFD,CPS  ,D,  ANS  , RE 
C 1 = C DS - 

$ = E F D 
FNO  = AN S 
AR  FA  = A LD  HA 

p7=3. 1 41 59 

S 2=  S / 2 . 

OMS  2=1.-$? 

A-S-P  = AUsl&  OMS? — — — — — — 

HlsClMP  1/2  . )<•  S?:*o^(  S/3  ) *FNO«-M  l.-S/B) 
H2  = OMS2  VH1 

DA=(  AU-AU  /FLOAT(  NN-1  ) 
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c 

c 

•€- 


C 

e 


?K  = lu-DA^LOAT<K-l) 

AR<  = AK /AU 
DE  LT  = 1 .-ARK 

TK=  ( AS P- AKi -i* 0MS2)  /H2 

L=  NN- ( K— 1 ) 

G(U  =A  RF  A*"  A K;i  M-$2  )?  { T K/  RET  A)  *»■'  ( AR  FA-1)* 

1EXP(  -(  TK /BETA)  mRFA)  /(HI*  BETA  ) 

ie< 6 ft  

IF(G(L).LT.l .0E-5  0 ) 50  TO  20 
10  CONTINUE 
GO  TO  40 

-?^-G^T+N-UE — — 

no  30  1=1  ,L 
3 0 0(  1)  =0.0 
40  RETURN 

en-u 

FUNCTION  OM  EXP  ( X ) 

B=-  50.  0 
C=1  .0F-20 

1E9ABS  (X l.-LT  .C)  GO  TO -5 

IF  (ABS(X)  .IT.  0.01  ) GO  TO  10 
IF(X.LT.B)  GO  TO  30 
OMEX  P = 1.  -EX  P (X  ) 

RETURN , — - 

5 0MEXP  = -0. 99999  9*X 
RETURN 

10  ONE  XP  = -0. 1 66666*  X*X*X- 0.499 999*  X'  X-0. 999999 *X 

RETURN — — ■ - ■ ^ : — - 

30  OMEX  P = 1 • 

RETURN 

END 

S UB-R-OUT-f  NB-O-SFU+l-r^T^T-NUl-m — — — 


01  MENSION  Y<  1)  ,Z(  1) 


HT  = . 3333333*  H 

I F { ND I M-  5)  7 » 8,  1 


1 SUM l = v ( 2 ) +Y( 2 ) 

SUM1=SUM1*SUMI 

SUM  1=9  T*'  ( Y(  1 ) * SUM  1 *-Y(  3 ) ) 

— AJX1  = Y(4  ) t-Y -4 4-1 — 

AUX1 =AUX 1*AUX1 

AUX1  = SUM1MT«M  Y(3)  + AJXUY(5)  ) 

AUX2=HT* ( Y ( 1 )*3.8  75MY(2  )+Y(  5)  ) *2 . 625*(  Y{  3)  * Y(  4)  ) * Y(  6)  ) 

— S U M-2-  = Y45L*Y  4 6) 

SUM2  = SUM2*SUM2 

S'JM2  = AUX2-9T*' (Y(4)+SJM2*Y(6  ) ) 

Z(  11=0. 

A U X = Y ( 3H-Y(  3 ) — 

AUX= AUX+  AUX 

Z(  2 ) = SUM  2-9  T1*  { Y(  2 ) *AUX  *Y(  4)  ) 

Z ( 3 ) = S UM  1 

7(4)  =SUM2 

IF(  NOT  M-6) 5,  5,  2 

INTEGRATION  LOOP 
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on 


2 00  4 I =7 ,ND I M,R 

sum  i = aux  1 

$UM2  = AUX  2 
AUX1=Y(I-1)  fY(  1-1  1 

— — AUX  1 = A UX  l +-AUX1 — — — 

AUX  1 = S UM 1 H T“~  ( Y ( I-2I+AUX1+Y  (I)  ) 

Z(I-2)  =SUM1 
IF  ( 1-NOIIM  13,6,6 

3-  AUX  2 =Y4  1 4+Y  HI 

AUX2  = AUX2*AUX2 

AJX2=SUM2HTHY(T-1)  + AUX2*Y(  TH  >1 

4 Z(I-l) =SUM2 

S—Z  (A}-!>FM — LHAUX4- 

Z ( M D I M ) = AUX  2 
RETURN 

5 Z(  X DIM-1  ) = SUM2 

I ( NO  I M 1 = AUX  1 — - 

RETURN 

C ENO  OF  INTEGRATION  LOOP 

C 

7 iPimm-3)  12,11,3 — 

c 

C MO  I M IS  EQUAL  TO  4 OR  5 

8 SUM2=l.l25f  HT*  ( Y(  IUYI  2UYI2)  fY(2l  U(  3)»V(3I  U(3UY(4)  ) 

S JM 1 * Y (2  M-Y424 — 

SUM1  = SUNHSUM1 

SU^I  = HT"(Y(  1 H-SUMlfY(  3 ) ) 

zm=o. 

AOXHYf  3 HYI  31 

AUX 1 = A UX 1 f A UX1 

Z ( 2 ) =SUM2-UTHY(?  ) *-AJX  1U  (4  11 
I F ( NDIM-5) 10,9,9 

9 AJXHY  444+Y444 

AUXl=AUXlf  A UX1 

Z(5  ) = SUM1*HT*(  Y(3  ) fAUX  UY(  511 
13  Z ( 3 ) =S  UM 1 

?f  44^-SUM? — — — 

R FT  JRN 


NO  I M IS  FOJAL  TO  3 

44  - S JM1  = UT-'*  (1 . 2 5^X14  1 +Y I 2 >+Y(  2 1-  .26-Y(3)  1 
SUM  2=  Y ( 2 I Y ( 2 J 
$UM2=SUM2*,SUM2 
7(3  )=HTMY<1  M-SUM2  4-Y(3  1 ) 

Z ( 4 4 *0-*- - - - - - — — 

7 ( 2 ) = SUM  1 


12  RETURN 
ENO 

F4X&T  I ON  FS4  A4 — 

C0MMON/POI/  CAi , NC , NR 
IF(A.LT.CAl)  F C=  0 . 
IF(A.GE.CAl)  FC=1. 


RETURN 

FNO 

//GO. SYS  IN  OD  * 
12  3.  0 

-423.0 

n 
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